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A CRITICAL  REVIEW  OF  TOE  NUMERICAL  SOLUTION  OF  NAVIER-STOKES  EQUATIONS 

by 

Sin- I Gieng 
Princeton  University 

Department  of  Aerospace  and  Mechanical  Sciences 
I.  INTRODUCTION 

This  article  concerns  the  various  practical  problems  of  using  high 
speed  electronic  computers  to  obtain  approximate  solutions  of  various 
fluid  flow  problems,  not  with  the  mathematical  techniques  of  solving 
the  Navier-Stokes  equations  through  difference  approximations  in  gener- 
ality. The  boundary  conditions  are  as  important  as  the  partial  differ- 
ential equations  in  the  mathematical  formulation  defining  a given  physi- 
cal problem.  There  are  conplicated  practical  problems  of  discretizing 
the  differential  formulation  (of  both  the  differential  equations  end 
the  boundary  conditions)  into  appropriate  difference  formulation  for  its 
numerical  solution.  There  are  quite  a few  unusual  aspects  in  such  at- 
tempts. 

Fluid  dynesdeists  usually  ignore  the  question  of  convergence  in  the 
asymptotic  differential  approximations  through  perturbation  arguments. 
They  often  consider  the  computational  solution  of  the  resulting  system 
of  ordinary  differential  equations,  as  routine,  even  though  it  may  be 
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very  tedious,  and  often  ill  posed  especially  for  the  multi -eigen -value 


Under  difficult  circumstances,  heurisitc  local  treatments 


are  often  introduced  to  reredy  the  situation,  or  a slightly  different 
original  approach  will  be  esqployed.  Such  a hit  and  miss  approach  has 
been  carried  over  to  the  direct  concur ational  solution  of  the  partial 
differential  equations  system  of  fluid  dynamics.  They  lead  to  much 
□ore  disappointment  and  often  tiros  with  more  serious  consequences. 
While  it  may  not  be  crucial  to  appreciate  the  mathematical  details, 
it  is  important  to  be  aware  of  the  ieplicaticns  of  some  fundamental 
mathematical  results  concerning  the  difference  approxi nations  of  a 
partial  differential  equation.  Accordingly,  a brief  review  of  these 
mathematical  aspects  will  be  outlined  prior  to  the  discussion  of  the 
practical  art  of  numerically  integrating  the  partial  differential 
equations  system  of  fluid  dynamics. 

Within  the  continuum  description,  the  fluid  will  be  considered  to 
be  homogeneous  and  to  possess  two  independent  thermodynamic,  or  state 
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variables,  i.e.,  the  density  p and  the  internal  energy  e per  unit  mass. 
There  is  an  algebraic  equation  cf  state,  relating  the  thermodynamic 
pressure  p to  its  density  p end  internal  energy  e as  p - p(p,e).  Let 
ui  be  the  velocity  vector  of  a fluid  element  with  i ■ 1,2  end  3 in  a 
three-dimensional  space  x^ . u j,  p,  and  e are  the  five  dependent  variables 
and  will  be  considered  as  functions  of  x.^  and  t.  The  Eulorian  descrip- 
tion of  the  change  of  these  dependent  variables,  is  the  set  of  five 
partial  differential  equations,  expressing  the  conservation  of  mass 
momentum  and  energy  written  here  in  divergence  form  as: 

^*^[puj]"°  a. 


3(pUi)  a r . 1 

Tt  - * 75$  [fVj  * p5ij  - Ti  j J * 0 

[puj(e*Vi/2)  ♦ puj'qj'uiTij]  • 0 


When  the  surface  stress  is  related  linearly  to  the  strain  rate  as 

3ui  Sui  3 *». 

Tij  * ^ J " 2 U)  3S7 

and  when  the  heat  transfer  vector  q^  is  linearly  related  to  the  temper- 
ature gradient  as 


qj  ‘ ’ k 


IsjT 


the  system  of  equaticns(i. 1-1.3) will  be  referred  to  as  the  Navier-Stokes 
equations  for  a coapressible  fluid.  The  Prandtl  nusber  Pr  and  the  specific 
heat  ratio  y are  properties  of  the  fluid  and  are  both  of  0(1).  The  shear 
viscosity  coefficient  p is  assumed  to  be  a known  algebraic  function  of 
temperature  (or  internal  enorgy).  The  bulk  viscosity  coefficient  is 
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often  taken  as  zero  or  otherwise  absorbed  in  u.  In  dimensionless  form 

a Reynolds  number  may  be  defined  in  terms  of  some  characteristic  length 

L and  velocity  U as  Re  « p U L t\i  where  subscript  o indicates  that 
o o oooo  o 

the  quantity  is  to  be  evaluated  at  some  reference  state.  For  most  fluid 
dynamics  applications,  the  Reynolds  number  is  very  large. 

The  divergence  form  of  the  Mavier-Stokes  equations  (1.1) -(1.3)  may 
be  written  as 


3v 

Ft 


1 3v 


(1.6) 


posed  at  an  initial  value  problem  for  the  vector  unknown  v,  having  the 
five  scaler  components  p,  pu^ , and  e.  F is  the  flux  of  v,  given  by  the 
nonlinear  quantities  in  the  square  brackets  of  Equations  (1.1) -(1.3). 

When  physically  meaningful,  initial  and  boundary  data  are  prescribed.  Equa- 
tion (1.6)  is  expected  to  siva  a satisfactory  description  of  the  tesporal 
development  of  the  flow  field  at  later  times.  This  expectation  is  mathe- 
matically justifiable.  The  integration  of  this  equation's  system  is  needed, 
for  example,  in  weather  forecasting  and  in  the  determination  of  the  temporal 
development  of  blast  waves,  hurricanes,  or  turbulent  fluctuations,  etc. 
(where  the  gravitational  field  and  the  coriolis  forces  are  included  where 
necessary).  In  most  aeronautical  applications,  the  steady  state  (or  the 
quasi-steady  state)  problems  are  more  often  of  primary  interest  where  the 
teatporal  dependence  is  neglected.  Thus  Equation  (1.6)  becomes 


^F(V.  -0 


3xk 


(1.7) 


which  is  to  be  solved  as  a boundary  value  problem.  The  boundary  con- 
ditions mu3t,  of  course,  be  independent  of  time.  But  it  is  not  clear 
how  such  boundary  conditions  should  be  specified  to  provide  the  required 
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steady  state  solution  or,  indeed,  any  solution  at  all.  Physical  intuition 
often  provides  meaningful  guidance  but  not  all  what  is  needed. 

The  stress  and  the  heat  conduction  terms  give  rise  to  the  second 
and  the  highest  order  partial  derivatives  with  coefficients  proportional 
to  Re”1.  This  steady  state  Navier-Stokes  Equation  (1.7)  generally  assumes 
an  elliptical  behavior.  When  Re  becomes  large,  the  flow  field  may  be 
divided  into  sii>- regions.  In  the  region  sufficiently  far  away  from  any 
solid  boundary,  the  inviscid  approximation,  obtained  by  dropping  terms 
in  (1.6)  or  (1.7)  containing  Re'1,  is  a valid  approximation,  known  as 
tho  Euler's  Equation 


3y 


1 1 * 3i7  F(v)  " 0 

j 


PM  - o 


(1.8a) 
(I. 8b) 


The  time  dependent  inviscid  Equation  (1.8a)  remains  hyperbolic  end  is 
posed  as  an  initial  value  problem  as  is  Equation  (1.6).  The  steady  state 
Equation  (1.8b),  however,  can  be  purely  elliptic  (subsonic)  or  purely 
hyperbolic  (supersonic)  or  mixed,  (i.e.,  with  both  elliptic  and  hyperbolic 
regions,  the  boundary  of  which  will  depend  on  the  solution  and  are  not 
known  before  hand  such  as  in  the  supercritical  transonic  inviscid  flow 
problem) . 

In  the  regions  near  the  solid  boundary,  or  near  where  there  is  a 

large  shear  stress  or  heat  conduction,  some  or  all  of  the  stress  terms 
X 3v 

contained  in  F(v,  ) have  to  be  kept  with  the  other  terms  despite 

the  large  Re.  If  this  viscous  region  should  extend  along  a coordinate 
surface  (Xj.Xj)  such  that  the  lateral  extent  (along  x^)  of  this  viscous 
layer  is  small  compared  with  its  physical  extent  along  the  (x1,x2)  sur- 
face, then  Prandtl's  boundary  layer  theory  applies.  Only  the  highest 


— i.— ■— 1 - - 
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order  partial  derivative  in  this  lateral  direction  (32/'dx*)  will  sur- 
vive the  limit  of  very  large  Re.  This  asymptotic  limit  at  large  Re  gives 
the  bowdary  layer  equations  which  are  parabolic.  However,  not  all 
viscous  layers  are  sheet  like  and  can  be  treated  by  Prandtl's  boundary 
layer  approximation.  For  such  viscous  layers  like  the  near  wake  and  the 
interaction  region  between  a shock  wave  incident  on  a boundary  layer,  the 
full  Navier  Stokes  Equation  (L7)  will  have  to  be  used  and  the  problem  be- 
comes elliptic  at  least  in  a significant  portion  of  the  flew  field  of 
interest. 

The  change  of  the  mathematical  character  of  the  flow  field  in  differ- 
ent regions  when  the  Reynolds  number  is  large  is  both  a blessing  and  a 
cause  of  concern.  It  is  a blessing  that  enabled  the  development  of  fluid 
dynamics,  historically,  in  the  form  of  the  inviscid  or  perfect  fluid 
theory  and  of  the  boundary  layer  theory.  But  it  is  also  the  fundamental 
difficulty  in  the  analysis  of  the  mixed  flow  regions,  characteristic  of 
mo3t  of  interaction  flow  problems.  Now  there  are  significant  differences 
in  the  numerical  integration  of  the  three  types  of  partial  differential 
equations.  The  method  that  has  been  proved  to  be  successful  for  one  type 
need  not  be  so  for  the  other.  It  is  therefore  iiqportant  to  recognize  the 
type  of  the  partial  differential  equation  before  formulating  a difference 
approximation  for  its  numerical  integration.  Clearly  then  there  are 
difficulties  in  the  numerical  integration  of  mixed  problems.  Such  diffi- 
culties are  quite  different  from  those  encountered  in  the  asymptotic 
analysis  of  mixed  flow  problems.  In  a few  examples,  they  have  been 
successfully  resolved  with  appropriate  cautionary  measures.  But  there 
is  no  theorem  to  guarantee  its  success  in  other  problems. 


> 


It  sight  be  asked  that  if  the  elliptic  steady  state  Navier-Stokes 
Equations  Q . 7)  could  be  integrated  for  a given  finite  but  large  Re,  why 
should  the  difficulties  arising  from  the  asymptotic  limit  of  Re  » con- 
cern us.  An  obvious  answer  nay  be  that  the  asymptotic  fora  of  the  partial 
differential  equations  system  is  such  simpler.  A more  fundamental  reason 
is  that,  at  finite  but  large  Reynolds  numbers,  the  asymptotic  behaviors 
of  the  flow  in  different  regions  bear  strongly  on  the  appropriateness  of 
the  difference  formulation  ar.d  the  numerical  integration  cf  the  Navier- 
Stokes  equations  when  the  resolution  (or  the  nimher  of  meshes  per  linear 
dimension  of  the  field  of  computation)  is  Sv..*rely  limited. 

For  a flow  problem  in  throe  space  dimensions,  an  average  of  30 
meshes  per  linear  dimension  will  give  rise  to  3 x 10*  nodel  points;  end 
will  need  1.5  x 10^  words  of  storage  spaces  for  the  5 unknowns  at  each 
point.  Such  storage  spaces  should  preferably  be  provided  in  the  core 
cf  the  computer  imit  for  ready  access.  Such  a requirement  will  stretch 
the  core  memory  capacity  cf  most  of  the  currently  available  large  com- 
puters such  as  CDC  6600  or  IEM  360-91.  The  solution  of  the  full  Navier 
Stokos  Equation  Ci.7)  for  a well-posed  boundary  value  problem  will  need 
hours  of  computation  in  such  machines.  The  parallel  computers  in  the 
stage  of  advanced  development  like  ILLIAC  IV  and  the  STAR,  can  neither 
promise  much  improvement,  in  this  regard.  To  extend  the  core  memory 
capacity,  a hierarchy  of  external  storages  will  be  providod.  Frequent 
reference  to  such  external  storages,  however,  will  greatly  increcse  the 
time  required  for  data  management  because  of  the  slow  input-output  devices 
connected  to  the  central  processing  unit.  This  problem  is  particularly 
aggravated  in  the  aforementioned  parallel  computers  where  the  promised 
large  gain  in  arithematic  speed  can  be  obtained  only  for  specific  inodes 
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of  "parallel"  or  "vector"  computations  in  which  a huge  amount  of  data 
mist  be  properly  processed  and  continuously  fed  into  the  ari thematic 
uait(s) . 

The  concept  of  parallel  use  of  an  array  of  aini-conputers  might 
appear  to  relieve  such  a difficulty  associated  with  the  specific  mode  of 
high  speed  arithematic  operations.  The  benefit  is  likely  to  be  illusory, 
however,  at  least  for  the  present  applications.  It  simply  transfers 
onto  the  users  the  tremendous  problem  of  optimal  coordination  of  the 
operations  of  the  array  of  mini -computers  and  the  problem  of  data  manage- 
ment among  the  di versed  "internal"  end  "external"  storege  facilities. 

The  u3ers  are  ill-equipped  with  the  expertise  of  the  computer  scientists 
who  designed  the  softwares  managing  the  businesses  of  the  central  pro- 
cessing unit  of  the  large  computers.  The  users  will  be  left  to  derive 
whatever  speed  advantage  each  individual  program  may  provide.  It  should 
bo  remeahered  that  without  the  order(s)  of  magnitude  increase  of  the 
"overall  processing  speed"  of  the  computer (s ) , the  increase  in  the  nua- 
bor  of  mesh  points  in  a linear  dimension  for  the  integration  of  a prob- 
lem can  easily  escalate  tho  computer  processing  time  from  hours  to  days 
and  months.  It  appears,  even  projected  slightly  into  tho  future,  that 
no  more  than  a couple  of  hundred  mesh  points  per  linear  dimension  can 
possibly  be  considered  in  the  integration  of  the  system  of  equations  of 
fluid  dynamics.  Kith  such  a limitation  on  the  resolution  of  corputational 
solutions,  the  integration  of  the  full  Havior-Stok.es  equation  for  the 
flow  field  over  a vehicle,  for  example,  (the  external  flow  problem)  seems 
futile.  The  field  of  computation  is  much  like  the  test  section  of  a 
win’dtunnel . Without  the  "full  scale"  facility,  coiqputers,  like  wind- 
tinnels,  should  be  used  at  present  to  study  only  the  components  or  the 


local  flow  fields.  For  such  purposes,  various  asymptotic  forms  of  the 
Navicr  Stokes  equation's  system  should  be  used  for  different  parts  of 
the  flow  field.  The  full  Ncvier-Stokes  equation's  system  will  be  called 
upon  only  for  the  study  of  those  flow  problems  that  cannot  be  consistently 
treated  with  the  simplified  flow  equations,  including  notably  those  mixed 
interaction  flow  problems.  Under  the  practical  limitations  of  resolution 
and  of  computer  time,  it  is  particularly  important  to  delineate  the  vary- 
ing nature  of  the  flow  regimes  in  the  different  treatments  of  the  Navier- 
Stokes  equations  end  to  consider  the  various  possibilities  of  how  the 
boundary  conditions  may  be  implemented.  To  quote  the  1S60  statement  of 
Forsythe  and  Wasow:  "The  numerical  solution  of  partial  differential 

equations  is  no  easy  matter.  Almost  every  problem  arising  out  of  the 
physical  sciences  requires  original  thought  arid  modifications  of  exist- 
ing methods".  This  statement  is  equally  true  today,  and  particularly  so 
for  the  type  of  flow  problems  under  consideration  here. 


II.  FUNDAMENTAL  CONCEPTS 


Consider  now  the  problem  of  solving  a partial  differential  equation 
subject  to  a set  of  initial  and  boundary  data  through  numerical  integration. 
A difference  formulation  is  obtained  by  replacing  the  differential  coeffic- 
ients by  appropriate  difference  quotients  as  an  approximation  to  the 
differential  problem.  There  will  be  some  "errors"  in  the  approximate 
formulation  both  in  the  equation  and  in  the  initial  and  boundary  data. 

When  these  "errors"  vanish  as  the  mesh  sizes  At  -*■  0 6 Ax  -*■  0 in  some  man- 
ner, the  difference  approximation  is  said  to  be  "consistent"  with  the 
differential  problem.  The  solutions  of  this  difference  formulation  pro- 
vide a sequence  of  approximate  solutions,  which,  in  the  limit  of  At, Ax  0, 
is  supposed  to  "converge"  to  the  solution  of  the  differential  problem  in 
some  sense;  i.e.,  the  "error"  of  the  solution,  as  a measure  of  the  depar- 
ture of  each  member  of  the  sequence  of  approximate  solutions  from  the  solu- 
tion of  the  differential  problem,  tends  to  "zero".  This  convergence  is, 
however,  not  guaranteed  for  a consistent  approximate  difference  formula- 
tion. Various  aspects  will  be  considered  in  the  following  sections. 

2.1  Well-Posed  Differential  Problem 

The  differential  problem  should  not  only  possess  a unique  solution, 

•» 

but  also  possess  "neighboring  solutions",  whether  the  problem  is  to  be 
integrated  analytically  or  numerically.  This  means  that  when  the  initial- 
boundary  data  is  slightly  perturbed,  the  differential  problem  should  still 
provide  a solution,  which  hopefully,  departs  from  the  unperturbed  solution 
of  tho  problem  only  slightly.  This  is  primarily  a physical  requirement 
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expounded  by  Hadamard  if  a mathematical  problem  is  to  describe  a physical 
situation  appropriately.  Mathematically  speaking  the  solution  of  the 
differential  problem  is  said  to  vary  continuously  with  the  data;  and  the 
differential  problem  is  said  to  be  "well  posed".  Among  other  things,  a 
given  partial  differential  equation  is  well  posed  only  when  the  boundary 
conditions  are  properly  specified.  For  example,  the  Laplace  equation, 
in  two  variables  x and  y. 


3Ju  3*u 

♦ yj?r 


0 


(2.1) 


is  well  posed  when  the  values  of  the  function  u(x,y)  is  specified  on  a 
closed  boundary  enclosing  the  domain  of  x.y  of  interest,  (Dirichlet  Problem). 
Now  the  function 

u(x,y)  ■ n"a  sin  nx  cosh  ny  (2.2) 


is  an  exact  solution  of  the  Laplace  ec.uation  with  the  initial  data 


u(x,0)  * n a sin  nx 
ly  (x,0)  - 0 

This  set  of  initial  data  is  small  everywhere  on  x with  a > 0 and  n suffic- 
iently large.  Now  if  the  Laplace  equation  is  to  be  solved  when  u(x,0)  and 

•§y  (x,0)  are  specified,  then  a small  perturbation  of  the  initial  data  can 
introduce  the  perturbBticns  of  the  type  (2.2)  onto  the  solution  of  the 
problem.  Thi3  perturbation  (2.2)  is  not  small  in  the  immediate  neighborhood 
of  y ® 0 despite  the  small  error  of  the  initial  data  when  n is  large  and 
a is  positive.  While  the  perturbation  (2.2)  does  vanish  at  y = 0 for  any 


value  of  n including  n ■*  °®,  the  value  of  u(x,y)  given  by  (2.2)  at  some 
small  but  finite  value  of  y,  bocomcs  infinitely  large  as  n + Thus 
the  Laplace  equation  is  not  "well  posed"  or  "ill-posed"  when  u(x,0)  and 

(x,0)  are  specified,  (Cauchy  Problem).  If  we  should  proceed  to  inte- 
grate this  "ill-posed"  problem,  the  perturbed  initial  data  is  expected  to 
contain  components  like  (?.2i  and  the  numerical  solution  will  r.ot  converge 
to  the  desired  solution  even  if  Ax  -*■  0 (i.e.,  n «•>) . 

If  the  gradient  of  u(x,y)  is  specified  over  a closed  boundary  (Neumann 
Problem)  or  if  either  the  gradient  or  the  value  of  u(x,y)  is  specified  over 
a closed  boundary,  the  problem  of  solving  the  Laplace  equation  is  well  posed 
if  some  integral  conditions  are  met.  Ill-posed  problems  will  result  other- 
wise, i.e.,  either  when  the  Dirichlet  or  Neumann  conditions  are  specified 
only  on  an  open  boundary,  or  when  the  Cauchy  condition  is  used  anywhere. 

Ihis  statement  is  applicable  to  elliptical  partial  differential  equations 
in  general.  The  parabolic  equations  sre  well  posed  under  similar  conditions 
but  only  on  an  "open"  boundary  and  when  integrated  in  the  "positive"  direction. 
The  hyperbolic  problems  are  well  posed  only  when  the  Cauchy  conditions  are 
specified  on  appropriate  "open"  portion  of  the  boundary.  It  becomes  then 
difficult  to  specify  the  boundary  conditions  that  will  render  a mixed 
differential  problem  well  posed  before  attempting  to  formulate  a difference 
approximation  of  the  problem  for  numerical  integration.  From  this  point  of 
view,  the  algebraic  convexities  of  the  full  Navier-Stokes  equations  system, 
either  for  the  time  dependent  hyperbolic  problem  (1.6)  or  for  the  steady 
state  elliptical  problem  (1.7),  may  well  be  tolerated  to  facilitate  the 
formulation  of  a well-posed  problem. 


2.2  We  11 -Posed  Differer.ee  Problem 
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It  is  not  only  that  the  differential  problem  should  be  well  posed 
for  a specific  or  a selected  class  of  initial  data,  but  also  that  the 
difference  problem  should  be  well  pcscd  snd,  for  a more  general  class  of 
initial  data,  to  provide  a convergent  numerical  solution.  This  is  because 
the  perturbations  implicit  in  the  numerical  solution  of  the  approximate 
difference  formulation  need  not  fall  within  the  class  of  the  initial  data 
for  which  the  differential  problem  is  well  posed.  The  function 

u(x,t)  ■ exp  [ia(x  ♦ t)]  (2.3) 

satisfies  the  first  order  hyperbolic  equation 

Jt  u ' ‘55T  u = 0 (2,4) 

with  the  initial  value 

u(x,0)  = exp(iax). 

The  complex  notation  with  i =\^T"is  used  here  for  simplicity,  to  mean 
that  both  the  real  and  the  imaginary  parts  of  the  expressions  should  be 
valid  simultaneously.  A wide  class  of  function  u(x,0)  can  be  formed  by 
superposing  various  trigonometric  initial  data  corresponding  to  various 
choices  of  the  values  of  the  constant  a.  Each  component  can  possess  an 
arbitrarily  assigned  amplitude.  By  summing  up  the  component  solutions  of 
different  o’s,  the  solution  of  the  problem  with  the  generalized  initial 
data  is  obtained.  Any  numbers  of  the  component  solutions  can  be  perturbed 
with  a correspondingly  small  perturbation  on  the  solution.  The  differential 
problem  is  thus  well  posed. 
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Stqppose  now  the  forward  tiae  centered  space  difference  algorithm  is 
used  to  provide  a difference  approximation  to  (2.4)  as: 


ur1  - \f!  u"  , - u" 


1L-—L1 

2Ax 


U°  ■ exp  (ictx) 

where  ■ U(jAx,nAt). 

An  exact  solution  of  the  difference  problem  (2.5)  is 

U*?  « U(jAx,nAt)  = (1  ♦ i jj-  sin  aAx)n  exp(iox)  (2 

where  n ■ t/At.  In  the  limit  of  At  -*•  0 S Ax  -►  0,  the  difference  solution 
U^,  given  as  (2.6),  converges  uniformly  to  the  solution  u(x,t),  given 
as  (2.3),  for  the  differential  problem  (2.4).  The  same  holds 
true  for  all  the  components  and  for  their  sum  with  a generalized  initial 
data.  Now  when  the  difference  problem  (2.5)  is  computed  for  any  small 
At  and  Ax,  the  computation  is  always  unstable  (as  is  well  known).  It  is 
apparent  that  some  components  of  the  perturbations  introduced  by  the 
computation  of  the  difference  form  (2.5)  cannot  be  represented  by  the 
trigonometric  data  and  grow  out  of  bounds. 

The  Euler's  equation  (1.8a)  for  inviscid  gas  dynamics  is  easily 
cast  into  Cauchy-Kowaleski  type  quasi-linear  hyperbolic  equations  system 


. AC)  §H- . 0.  I 

3F 

where  A(u)  » . If  the  initial  data  of  u(x,t  « 0)  ■ f(x)  and  A(u)  are 

analytic,  then  the  solution  u(x,t)  for  all  x and  t is  analytic.  The 
requirement  of  the  analyticity  of  the  initial  data  might  not  appear  to 
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to  be  very  restrictive  in  view  of  the  Weierstrass  approximation  theorem. 

True  that  any  continuous  function  within  a closed  interval  can  be  approxi- 
mated  arbitrarily  closely  by  enalytic  functions , including  polynomials 
and  sinusoidal  functions.  But  arbitrarily  close  approximation  of  the 
initial  data  does  not  promise  the  arbitrarily  close  approximation  of 
the  solution.  The  two  examples  (2.1)  and  (2.5)  given  above  illustrate 
this  point  both  for  the  differential  and  the  difference  equations. 

Equation  (1.8a)  or  (2.7)  for  inviscid  gas  dynamics  is  a well-posed 
problem  for  a fairly  broad  class  of  initial  data.  Even  if  it  is  presumed 
that  a consistent  difference  approximation  possesses  a solution  that  con- 
verges uniformly  to  the  solutions  of  the  differential  problem,  stable 
computation  is  not  guaranteed.  The  instability  of  the  computation  is 
attributed  to  the  fact  that  the  perturbations  on  the  initial  data  intro- 
duced by  the  computational  procedure  is  beyond  the  class  of  perturbations 
ejq>ressible  in  terms  of  the  piecewise  enalytic  data.  For  a difference 
problem  to  be  well  posed,  its  solution  oust  be  continuous  with  a much 
wider  class  of  perturbations  on  the  initial  data.  This  is  the  crux  of 
the  concept  of  computational  stability. 

Computational  stability  in  general  calls  for  the  boisidcdness  of  all 
the  perturbations  in  the  computed  solution.  Then  when  the  magnitudes  of 
the  perturbations  in  the  initial  data  are  made  arbitrarily  small  in  the 
limit  of  vanishingly  small  mesh  sizes,  the  resulting  perturbations  in 
tlie  computed  solution  will  likewise  vanish.  The  computed  neighboring 
solutions  based  on  a consistent  difference  formulation  will  then  converge 
to  the  solution  of  the  differential  problem;  i.e.,  stability  and  consistency 
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means  convergence.  This  is  the  essence  of  the  equivalence  theorem  of 
Lax.  The  success  in  obtaining  a convergent  approximate  solution  through 
the  computation  based  on  a given  difference  formulation  therefore  lies  in: 


(i)  the  consistency  of  the  difference  formulation  with  the 
well-posed  differential  problem, 

(ii)  the  stability  of  the  difference  formulation. 

Kore,  the  difference  formulation  means  collectively  all  the  difference 
relations  connecting  the  values  of  functions  at  different  time  levels  and 
at  all  the  mesh  points  in  the  interior  of  and  on  the  boundary  of  the 
field  of  computation. 


2.3  Computational  Stability 

Computational  stability  is  a characteristic  of  a set  of  difference 
equations,  not  of  a difference  algorithm  how  a differential  coefficient  in 
the  differential  equation  is  to  be  replaced  by  a difference  quotient. 

It  is  incorrect  to  refer  to  an  algorithm  as  stable  or  unstable.  The  same 
algorithm  when  applied  to  different  differential  equations  can  lead  to 
different  difference  equations  with  entirely  different  stability  character- 
istics. Thus,  the  forward  time  and  centered  space  difference  algorithm 
when  applied  to  the  simple  wave  equation  (2.4)  leads  to  an  always  unstable 
difference  equation  (2.5).  Mien  the  same  algorithm  is  applied  to  the 
heat  diffusion  equation,  the  resulting  difference  equation  is  stable  if 

s « < i.  , Some  simple  exanples  are  given  in  Tables  I and  II. 

fix'  — 2 

Slightly  different  algorithms . applied  to  the  same  differential  equa- 


tion  may  yield  difference  equations  with  quite  different  stability  be- 
havior. For  the  simple  wave  equation  (2.4),  the  forward  time  and  the 
backward  space  difference  algorithm  will  yield  an  always  unstable  difference 


equations 


U1?*1  - if  if  - if.  , 

1 Tt  2 •cJKJ~-  • 0 

The  forward  tine  and  the  forward  spatial  difference  algorithm  will  pro- 
vide the  difference  equation 

u P*1  - iP  , - iP 

-c  1 Vx  ■ 0 <2-9) 


that  will  be  stable  if  C At/Ax  £ 1.  And,  as  mentioned  previously,  the 
forward  time  and  centered  space  difference  algorithm,  Equation  (2.5),  is 
always  unstable.  The  choice  of  difference  algorithm  for  discretization  to 
obtain  a stable  difference  equation  is  not  trivial. 

For  a slightly  more  complicated  equation,  the  situation  is  considerably 
more  complex.  A partial  differential  equation  of  higher  order  may  be 
written  as  an  equivalent  system  of  lower  order  partial  differential  equa- 
tions. (Contrary  to  the  situation  of  higher  order  ordinary’  differential 
equations,  this  is  not  generally  true  for  partial  differential  equations). 

A partial  differential  equation  representing  a physical  principle  may  be 
written  in  different  but  equivalent  forms  in  terms  of  different  subsidiary 
variables.  When  the  same  difference  algorithm  is  applied  to  discretize 
these  equivalent  differential  forms,  the  resulting  difference  equations 
are  not  equivalent  and  may  possess  widely  different  stability  behavior  in 
coaqmtation.  Consider  the  simplest  case  of  the  second  order  wave  equation 


(2.10) 


which  is  equivalent  differentially  to  a system  of  two  first  order  wave 
equations.  We  may  write  the  system  in  terms  of  different  variables  as: 
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When  forward  time  and  centered  space  difference  algorithm  is  used,  the 
following  difference  equations  system  results 


/ *n+1  - 0n 

vr ; „ *J-r2J  * “j'1 


(2.12) 


and 


/ V?*1  - V1?  U1?  , - U?  . 
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I U0*1  - u”  V11  - V11 


(2.13) 


2Ax 


The  system  (2.12)  is  always  unstable  for  any  choices  of  At  end  Ax  (easily 
verified  by  v.  Neumann  Analysis)  while  the  system  (2.13)  is  stable  if 
C At/Ax  <.  1.  Note  also  that  the  similar  difference  equation  (2.5)  for  the 
first  order  wave  equation  (2.4)  is  always  unstable.  Thus,  it  is  not  a 
matter  of  trivial  consequence  to  rewrite  a partial  differential  equation 
into  equivalent  but  different  forms  before  discretization  with  the  same 
difference  algorithm. 

The  equations  of  fluid  dynamics  represent  the  three  conservation  laws 
of  mass,  momentum  and  energy.  They  can  be  expressed  in  terms  of 
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a great  number  of  dependent  and  independent  variables  or  in  terms  of 
particular  combinations  of  such  variables  and  in  different  coordinate 
systems.  The  second  order  equations  may  be  split  into  first  order  sys- 
tems. (For  the  moment,  the  question  of  nonlinearity  is  put  aside.)  For 
all  these  varied  forms  of  equivalent  systems  of  partial  differential  equa- 
tions (equivalent  in  the  sense  of  physics  and  mathematics)  a given  differ- 
ence algorithm  will  give  correspondingly  varied  difference  forms  with 
quite  different  stability  and  other  computational  behaviors. 

The  complete  difference  formulation  of  a fluid  dynamics  problem 
will  call  for  the  discretization  of  the  boundary  condition  in  addition  to 
the  set  of  differential  equations.  The  set  of  difference  relations  connec 
ting  the  values  of  various  functions  at  mash  points  neighboring  the  boun- 
dary is  generally  different  from  the  set  cf  recursive  relations  for  the 
interior  points  derived  from  the  differential  equations.  This  boundary 
set  of  difference  relations  may  be  unstable  while  the  recursive  differ- 
ence relations  for  the  interior  points  are  stable.  Apparently  trivial 
modifications  of  the  difference  formulation  in  the  boundary  conditions 
thus  often  lead  to  substantial  changes  in  the  stability  behavior. 

In  view  of  such  a complicated  situation  and  of  the  frequent  experi- 
ence of  severe  computational  instability,  it  it  highly  desirable  to  be 
able  to  analyze  the  stability  behavior  of  a given  difference  formulation; 
but  there  is  no  simple  asans  available  except  the  so-called  "energy 
analysis".  The  "energy  analysis"  attempts  to  establish  a finite  bound  of 
the  solution  (over  the  entire  net  or  mesh  space)  in  some  suitably-chosen 
norm,  and  the  formulation  is  by  definition  stable.  When  such  a bound  is 
established,  the  proof  of  convergence,  existence,  and  uniqueness  follows 
trivially.  For  a nontrivial  boundary  value  problem,  such  a proof  is  very 


difficult  and  very  tedious  even  for  a simple  equation.  Such  proofs  are 
available  for  the  Navier-Stokes  equations  for  an  incompressible  fluid,  but 
only  for  the  periodic  boundary  conditions,  a case  which  is  really  not  that 
much  different  from  a pure  initial  value  problem.  With  rather  complicated 
boundary  conditions,  it  is  not  practical  if  not  impossible,  to  ascertain 
the  stability  property  of  a difference  formulation  of  a fluid  dynamics 
problem  via  such  an  approach.  At  present,  it  is  a practical  art  to  draw 
from  experiences  with  similar  problems  and  inferences  of  model  analysis 
in  formulating  the  recursive  difference  relations  for  the  interior  points. 

The  formulation  of  the  boundary  conditions  is  approached  on  an  individual 
basis  and  modified  where  necessary.  The  entire  program  is  then  tested  in 
actual  computer  computation  for  its  stability.  Considerable  work  will  be 
involved  before  stable  computation  is  achieved.  By  then  quite  a few  modi- 
fications may  have  been  introduced.  It  is  opportune  to  check  if  the  final 
difference  formulation  is  consistent  with  the  differential  problem  to  be 
solved  both  as  to  the  differential  equations  and  the  boundary  conditions. 

It  may  well#be  that  the  physical  boundary  conditions  that  should  be 
applied  are  quite  different  from  those  consistent  with  the  difference 
formulation  or  that  some  spurious  terms  might  have  been  introduced  into 
the  differential  equations  that  fail  to  vanish  in  the  limit  of  At, Ax  0. 

The  v.  Ncir*  >«r  stability  analysis  tor  the  local  linearized  model  will  most 
likely  iqpose  some  restrictions  on  At  and  Ax  for  the  computation  to  remain 
stable.  This  restriction  should  be  observed  by  all  the  approximate  solu- 
tions, as  successive  members  in  the  Cauchy  sequence,  converging  toward 
the  solution  of  the  differential  problem.  The  limit  process  in  the  t-x 
space  is  not  to  be  taken  in  any  arbitrary  manner.  This  restriction  should 
be  considered  while  investigating  the  consistency  of  the  difference  formulation. 


23  - 


Certain  difference  algorithms  are  often  referred  to  as  "unconditionally 
stable".  What  it  means  is  that  when  such  an  algorithm  is  used  to  discretize 
a certain  type  of  differential  equation  for  the  solution  of  pure  initial 
value  problems  or  of  periodic  boundary  value  problems,  there  will  not  be 
restrictive  conditions  on  the  choice  of  At  (or  of  the  iterative  steps  in 
the  solution  of  pure  boundary  value  problems)  for  a given  set  of  Ax,  accord- 
ing to  the  v.  Neuman  stability  analysis  of  the  linear  equations.  When  such 
an  algorithm  is  used  in  the  numerical  solution  of  non-periodic  boundary 
value  problems,  even  for  that  particular  type  of  equations,  computational 
instability  will  often  result  especially  for  complicated  boundary  condi- 
tions and  for  non-linear  equations.  Even  if  no  question  of  stability 
arise,  the  apparent  advantage  of  permitting  the  use  of  too  large  time  steps 
At  need  not  lessen  the  computing  time  while  inevitably  decreases  the  accur- 
acy of  the  computed  solution.  Indeed  it  is  advisable  under  the  circum- 
stances to  verify  the  consistency  conditions  for  both  the  equation  and  the 
boundary  condition. 

A case  to  illustrate  the  point  is  the  following.  The  integration  of 
the  simple  heat  diffusion  equation 


3u  32u 


with  the  formally  second  order  accurate,  centered  time,  centered  space 
algorithm  of  DuFort- Franks 1: 
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is  "unconditionally  stable"  for  any  positive  values  of  s * so  that 

At  can  be  made  as  large  as  Ax  or  larger  without  leading  to  computational 


instability.  Most  other  explicit  difference  algorithms  when  applied  to 
the  heat  diffusion  equation  will  impose  a stability  limit  like  s <_  V 2 . 

The  restriction  on  At  is  particularly  severe  at  small  Ax.  Now,  in  inte- 
grating Equation  (2.15),  it  is  tempting  to  use  as  large  a At  as  is  practi- 
cal, usually  comparable  to  Ax,  to  save  computing  time.  Indeed,  this  is 
often  credited  as  the  "merit"  of  the  Dufort-Frankel  scheme.  Equation 
(2.15)  is,  however,  consistent  with  the  heat  diffusion  equation  only  when 
0 as  Ax  •*  0.  Otherwise,  it  is  consistent  with  the  wave  equation, 
having  the  wave  speed  . 
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(2.15a) 


With  2^  * 0(1) » the  computed  solution  is  expected  to  display  waves  compar- 
able to  the  temporal  and  the  spatial  variations  of  the  solution,  and  there- 
fore loses  much  of  its  value  as  an  approximation'  to  the  solution  of  the 
diffusion  problem  (2.14).  Even  with  At/Ax2  = s <_  1/2  for  example,  the  com- 
puted solutions  will  still  display  oscillations,  albeit  at  smaller  ampli- 
tudes. The  mean  solution  (taken  over  the  waves)  is  neither  a meaningful 
approximation  of  the  solution  to  the  diffusion  problem  with  Dirichlet 
boundary  conditions.  If  a Neuman  boundary  condition  should  be  imposed, 
instability  will  result.  The  qualitative  statements  mentioned  here  should 
not  be  generalized.  The  simple  example  is  given  above  only  to  bring  home 
the  point  that  every  individual  problem  should  be  carefully  examined 
according  to  the  fundamental  principle.  Our  current  understanding  of  the 
numerical  integration  of  the  partial  differential  equations  does  not  warrant 
any  simple  generalization  to  be  applied  to  the  complicated  situations  of 
fluid  dynamics. 


III.  STABILITY  ANALYSIS 
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The  numerical  integration  of  the  Navier-Stokes  Equations,  as  an 
outstanding  example  of  complicated  partial  differential  equations  system 
is  expected  to  encounter  quite  serious  practical  difficulties.  Such  diffi- 
culties fall  into  three  basic  categories  which  may  not  be  all  independent: 

■ (i)  Computational  Stability  All  disturbances  will  remain  bounded 
in  the  computation.  Otherwise,  the  value  of  some  quantity  will  eventually 
become  so  large  as  to  l*e  beyond  the  capability  of  any  computer  and  no 
results  would  be  obtained.  Hence,  this  is  often  referred  to  as  computa- 
bility. 

(ii)  Convergence  Rate  - The  solution  at  some  later  time  T or  at  the 
asymptotic  steady  state  should  be  obtained  with  a reasonable  amount  of 
computational  work,  i.e.,  the  number  of  time  or  iterative  steps  in  the 
solution  must  not  be  too  large  and  the  computational  work  for  each  step 
not  excessive  so  that  results  can  be  obtained  within  a reasonable  amount 
of  time  ( and  hence  .cost) . 

(iii)  Accuracy  - The  solution  eventually  obtained  must  in  some  sense 
approximate  the  physical  results  in  question  for  it  to  be  useful.  The 
criterion  for  its  being  an  adequate  approximation  is,  however,  subject 
to  judgement.  The  accuracy  criterion  imposes  limitations  on  the  fineness 
of  the  resolution,  both  temporally  and  spatially,  which  in  turn  sets  the 
requirement  on  the  convergence  rate. 

Conputational  stability  is  clearly  the  most  pressing  problem,  since 
it  is  the  first  one  to  be  encountered  in  an  attempt  to  get  any  solution. 
Much  work  has  been  devoted  to  this  question.  As  is  explained  in  the  pre- 
vious chapter,  its  fundamental  nature  is  essentially  understood  but  there 
are  quite  a few  subtle  aspects  in  its  implementation  even  for  the  simple 
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examples.  The  practical  question  of  achieving  a stable  computation  for 
the  complicated  system  of  Navier-Stokes  Equations  is  expected  to  be 
formidable.  The  various  heuristic  approaches  that  promise  to  guide  the 
formulation  of  a stable  difference  problem  will  be  reviewed  in  the  follow- 
ing chapter.  Generally  speaking,  with  some  hard  work,  stable  computation 
can  usual lly  be  achieved  as  may  be  verified  in  actual  computation.  It 
is  inyortant,  however,  at  this  juncture  to  bear  in  mind  that  the  convergence 
rate  and  the  accuracy  of  the  formulation  should  not  be  seriously  compromised 
in  an  all-out-effort  to  achieve  stability  of  the  computation.  The  objec- 
tive of  the  confutation  is  to  obtain  valid  approximations  to  a given 
physical  problem.  In  the  following  review,  it  is  therefore  intended  to 
bring  out  primarily  the  mathematical  assumptions  and  their  physical  inpli- 
cations  of  various  approaches  when  they  are  applied  to  the  solution  of 
different  types  of  practical  problems. 

3.1  v.  Neumann  Stability  Analysis 

A vector  unknown  function  U(t,x^)  of  dimension  p is  to  be  calculated 
over  mesh  spacings  Ax^,  Ax^,  for  successive  increments  of  At  from  the 
initial  values  of  U(t*0,x^)  based  on  a system  of  linear  difference  equa- 
tions . 

The  general  form  of  the  linear  difference  relation  may  be  that  some 
linear  combinations  of  the  values  of  the  function  Un+*  at  a group  of 
neighboring  mesh  points  are  given  by  some  ether  linear  combination  of  U11 
at  various  neighboring  points.  If  only  the  Un+*  evaluated  at  a single 
mesh  point  is  involved  in  the  difference  equations,  the  unknown  values  of 
U"+1  at  any  given  mesh  point  can  be  determined  without  reference  to  the 
advanced  values  of  l/1**  at  other  mesh  points.  Such  difference  equations 
are  explicit.  If  the  advanced  values  of  Un+*  at  more  than  one  mesh  point 
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are  involved,  a set  of  recursive  difference  relations  written  for  all 
the  mesh  points  would  have  to  be  solved  simultaneously  so  that  the  ad- 
vanced values  of  all  the  mesh  points  in  the  entire  field  of  computation 
will  be  obtained  at  the  same  time.  Such  difference  equations  are  implicit. 
Sometimes  it  is  preferred  to  solve  simultaneously  for  the  advanced  values 
at  special  groups  of  mesh  points  in  succession,  such  as  by  rows,  by 
columns,  or  by  diagonals,  by  blocks  or  by  bands.  Such  difference  equa- 
tions are  partially  implicit  and  partially  explicit  by  nature.  The  organi- 
zations of  the  special  group  may  change  from  one  to  the  next,  and  such 
different  groups  are  often  applied  in  alternate  sequence  or  in  some  special 
orders.  They  are  then  referred  to  as  alternating  direction  methods.  The 
specific  details  of  the  difference  algorithms  that  may  be  employed  to 
represent  a differential  problem  is  indeed  very  numerous. 

If  all  the  coefficients  of  the  difference  equations  are  constant,  and 
if  the  system  of  equations  are  to  be  solved  under  periodic  boundary  con- 
ditions (or  under  the  presumption  that  tire  boundary  is  so  far  away  as  to 
exert  no  influnece  on  the  solution,  i.e.,  the  pure  initial  value  problems), 
the  solution  of  the  system  of  equations  can  be  extended  periodically 
beyond  the  field  of  computation  with  both  u”  and  U°+*  represented  by 
Fourier  series.  The  linearity  of  the  difference  equations  system  permits 
the  treatment  of  each  Fourier  component  separately.  Thus  by  substituting 
U by  V(kj)  exp  (ik^xj  into  the  system  of  difference  equations  and  can- 
celling the  common  factor  in  each  equation,  an  equation 

= HoVn(k.)  (3.1) 


A 
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results.  Here  i is  the  complex  number  to  represent  the  sinusoidal  func- 
tions with  wave  numbers  k^,  k^,  k^  in  the  x^,  x x^  directions  respec- 
tively. V(k^)  is  the  amplitude  of  the  particular  wave  component  under 
consideration.  Each  of  the  Fourier  components  may  be  considered  either 
as  a part  of  the  proper  solution  U or  as  a small  perturbation  (or  error) 
superposed  on  the  solution  U.  and  Hq  are  the  matrix  operators  de- 
pending on  the  constant  coefficients  of  the  difference  equations  and  of 
At  and  Ax..  On  the  assumption  that  can  be  inverted,  then  Equation  (3.1) 
becomes 


v"+1(k.)  * G(At,Ax, ,k .)Vn(k .) 


‘j 


where 


r y 


.-i. 


(3.2) 


G(At,Ax.,k.)  - (Hj)  *Ho 


Equation  (3.2)  tells  the  evolution  of  each  Fourier  component  either  as  a 
part  of  the  solution  or  as  a perturbing  error.  Accordingly  G(At,Ax^,k^) 
is  called  the  amplification  matrix  of  the  system  of  difference  equations. 
The  condition  that  the  solution  U should  be  uniformly  bounded  requires 
each  and  every  component  to  be  so  bounded.  Since  the  norm  of 
nil  igimivoii.  u is  necessary  and  sufficient  that  ||g|P  be  so  bounded 
for  all  wave  components  and  for  all  n = T/At  where  T is  the  time  period 
for  which  U is  to  be  calculated  with  some  choice  of  small  but  positive  At. 
Now  |G|f\>  Rn(At,AXj,  kj)  where  R is  the  spectral  radius  of  G,  i.e.,  the 
largest  eigenvalue  of  G.  Hence  for  such  initial  periodic  boundary  value 
problems,  the  v.  Neumann's  condition  follows  that  all  the  eigenvalues  of 
the  anp  lification  matrix  G be  < l+0(At).  It  is  a necessary  condition 
for  confutation al  stability.  The  eigenvalues  of  the  amp lification  matrix 
are  often  more  conveniently  obtained  by  direct  substitution  of  un**  ■ Xun 
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into  the  difference  equation  to  obtain  the  determinant  which  vanishes 
when  X takes  up  the  eigenvalues.  The  v.Neumam condition  becomes  sufficient 
for  the  stability  of  the  stated  problem  when  the  matrix  G is  normal.  Both 
this  sufficiency  aspect  and  the  additional  term  0(ht)  are  without  much 
practical  significance  for  the  present  consideration  as  will  shortly  be- 
come clear.  The  important  points  to  recognize  from  the  above  are  the 
physical  indications  of  the  various  conditions  under  which  the  v. Neumann 
stability  analysis  is  formulated. 

3.2  Local  Linearization 

The  application  cf  the  v. Neumann  analysis  for  the  stability  of  the 
numerical  integration  of  a system  of  nonlinear  partial  differential  equa- 
tions such  as  (1.7)  calls  for  quite  a few  important  additional  approximations: 

(i)  The  nonlinear  difference  (or  differential)  equation  is  linearized 
by  considering  the  solution  as  the  sum  of  a small  perturbation  (or  variation) 
superposed  over  the  local  solution  of  the  problem.  By  substituting  the 
perturbed  solution  into  the  difference  equation  and  keeping  only  the  terms 
involving  first  power  of  the  perturbation,  the  result  is  the  equation  of 

the  first  variation  with  coefficients  depending  on  the  solution  of  the 
differential  problem,  such  coefficients  will  vary  with  x and  t. 

(ii)  The  coefficients,  will  be  assumed  to  be  slowly  varying  so  that 
these  coefficients  can  be  evaluated  by  the  constant  local  values  at 
various  mesh  points.  The  system  of  equations  of  the  first  variation  then 
becomes  linear  and  with  constant  coefficients  for  the  first  variation.  It 
is  supposed  to  apply  locally  at  each  mesh  point.  The  coefficients  (and 
hence  the  difference  relations)  vary  from  mesh  point  to  mesh  point. 

(iii)  The  stability  behavior  of  the  computation  at  each  mesh  point 
is  independent  of  its  neighbors  and  can  be  considered  as  the  stability 
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problem  very  far  away  from  the  boundary  (it  cannot  be  periodic).  Thus 
the  v. Neumann  stability  analysis  may  be  applied  locally  to  find  the  local 
stability  limit  on  At. 

(iv)  The  local  stability  limit  is  detsimincd  at  every  mesh  point 
in  the  interior  with  the  local  computed  value  U°(x)  rather  than  the  genuine 
solution  u(x,nAt).  The  most  restrictive  of  the  local  stability  limits 
over  all  the  interior  points  will  then  be  taken  as  the  stability  limit  on 
At  for  the  integration  of  the  difference  problem. 

Such  a local  linear  stability  analysis  to  be  applied  to  fluid  dynam- 
ics problems  is  probably  what  led  v. Neumann  to  develop  the  Fourier  method 
for  the  constant  coefficient  linear  difference  equations.  This  method  is 
still  the  most  valuable  practical  tool.  It  should  be  noted  that  slight 
difference  in  the  linearization  procedure  can  lead  to  slightly  different 
linearized  equations  of  the  first  variation.  They  will  give  slightly 
different  local  linearized  stability  criteria.  Consider  the  following 
examp 16: 

It  * IxT  {s*5a) 

• It  (5°4  It*  (3-3b) 


with  the  initial  condition 

uQ  «=  u(x,t=0)  * ¥[v(x0-x)] 

and  the  boundary  conditions 

u(o,t)  = YlvCyt+x^] 
u(L,t)  - *f[v(vt-L+xo)  ] 
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The  following  is  a 'solution  representing  a running  wave  with  constant 
wave  velocity  v 

u(x,t)  =*  y[v(vt-x+XQ)]  (3.4) 

where  the  function  4*  is  given  implicitly  as  the  inverse  of 

7 Cii-Uq)4  ♦ — ■ uo(u-uo)3  ♦ 15  u2  (u-uo)2  ♦ 20  U(>3  (u-uq) 

♦ 5 u4  £n(u-uQ)  = v(vt-x+xQ)  (3.5) 

It  is  shown  in  Fig.  1 with  a relatively  sharp  front  and  approximately  a 
quartic  curve  far  downstreams.  It  may  be  interesting  to  note  that  the 
Equation  (3.3)  stands  as  a heat  diffusion  equation  with  variable  diffus- 
ivity  5u4,  rather  than  an  equation  describing  the  steady  propagation  of 
a r.cndecaying  wave. 

Let  the  Equation  (3.3)  be  discretized  with  forward  time  and  centered 
space  difference  algorithm  with  the  spatial  derivative  evaluated  in  part 
0 at  the  advanced  time  step  snd  the  other  part  (1-0)  at  the  original  time 
step. 

°j+1  • ♦ C1-6)[62(U5)]”}  (3.6) 

where  the  second  order  spatial  difference  operator 

[«2C  )]"  - ( )”+1  - 2(  )“  ♦ ( )^_j  (3.7) 

The  parameter  0 can  be  chosen  at  convenience.  This  is  a nonlinear  equa- 
tion. The  following  linearized  approximation 

(U5) j*1- (u5) j - SOI4)*}  (U^1-  Uj) 


[ 


- 32  - 


gives  the  equation  of  first  variation  of  (3.6)  as 

- «$>  - [«*;.,  a#  - «$.,> 

- 2('J4}"  (u"*1-  l/j)  . CU4)".,  «£}  - l/j.,)] 


At 

ST 


(3.8) 


Note  that  the  equation  is  linear  in  the  inknown  (U^*1-  U*J)  if  the  values 
of  the  function  U at  all  the  spatial  mesh  points  at  the  time  level  n are 
known.  Otherwise,  the  equation  would  have  retained  its  nonlinear  fora. 
Alternatively,  it  is  also  appropriate  to  linearize  in  many  other  ways. 

A particularly  single  one  is  to  approximate 


and 


[«W>]y*1-  [dV)]j 

S(l,4)j  [(62U)"41-  (62U)"J 

[«V)]5  « 5(U4)^  (62U)” 


(3.9) 


Then  the  equation  of  the  first  variation  of  (3.6)  becomes 

CUj*1-  uj)  - sfj£  (U4)"  [o#  - uj„)  - 2CUJ*1-  uj)  . (Uj:‘  - 

■ s sir  [*,-  2^  ♦ “5-ij  (s-*°) 

This  last  equation  is  indeed  the  same  as  what  would  result  if  the  effec- 

4 

tive  diffusivity  5u  in  Equation  (3.3b)  should  be  treated  as  a constant 
before  and  during  discretization. 

With  all  and  taken  to  be  constant,  the  v. Neumann 
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stability  analysis  for  Equation  (3.8)  will  require  the  following  absolute 
value  to  be  less  than  unity  for  all  wave  number  k. 


1 - (1-6)5 
1 ♦ 8s 


< 1 


where  s is  the  cosqplex  expression 


(3.11) 


s 


5(U  Jj  S>7 


[* 


- (a+0)  cos  kAx  - i(a-0)  sin  kAx 
with  a - (U^/U^)4 

and  6 - (U^/U*?)4 


The  restriction  on  the  value  of  can  be  computed  for  all  k from 
(3.11)  for  different  values  of  a and  0 at  every  interior  mesh  point.  This 
is  a very  tedious  process.  The  v. Neumann  stability  analysis  for  Equation 
(3.10)  leads  to  the  same  relation  (3.11)  but  with  a*0»l.  This  provides 
an  explicit  limit  on  that  when  0 < 1/2, 

5(U  ’ SxT  - 2(1-20)  (3.12) 

and  there  will  be  no  limit  on  if  0 >.  1/2.  This  is  the  well-known 
results  for  the  simple  heat  equation. 

To  test  the  usefulness  of  the  local  linear  stability  criterion, 
coag>utations  were  carried  out  at  0 ■ 0.4  and  • 0.001  with  Equation 
(3.8).  The  parameters  v and  UQ  were  chosen  as  vj^-  » 0.075  and  5U4  = 0.005. 

The  lest  value  is  much  less  than  2. SO  os  is  required  by  Equation  (3.12)  for 
local  cocputational  stability.  As  the  computation  proceeds,  the  values 
of  U increases  with  t over  the  entire  field  of  computation.  According 
to  the  local  stability  criterion  (3.12),  we  would  expect  instability  to 
appear  in  the  form  of  rapidly  increasing  amplitudes  of  oscil'  • f.ions  when 
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and  where  the  values  of  (U1? -u  ) exceed  (500)  ^ ^ 4.7.  This  was  what 

J o 

happened  as  is  illustrated  in  Fig.  (1  ).  The  computed  points  lie  very 

close  to  the  analytical  solution  except  at  the  foot  of  each  wave  front 
where  the  solution  undergoes  a rapid  change  and  in  the  region  where 
U^/Uq  > 5 and  whore  the  confuted  solution  oscillates,  signalling  the  on- 
set of  computational  instability. 

It  is  remarkable  that  the  simple  local  criterion  deduced  for  the 
difference  Equation  (3.10)  provides  highly  satisfactory  guidance  for  the 
integration  of  Equation  (3.8)  although  o and  S generally  differ  from 
taiity.  Where  the  stability  boundary  of  Equation (3. 8)  with  a 8 0 + 1 lies 
within  the  stable  region  of  Equation  (3.10),  the  local  linear  stability 
limit  need  not  even  be  "necessary"  at  such  interior  points.  Such  regions 
are  likely  to  be  small,  hewevor,  if  the  model  equation  (3.10)  in  the 
above  example  is  appropriately  chosen.  It  is  clearly  not  sufficient 
since  the  v. Neumann  stability  condition  itself  is  not  and  since  the  in- 
fluence of  the  boundary  conditions  on  computational  stability  is  yet  to 
be  investigated.  Nevertheless,  the  local  linear  stability  analysis  does 
appear  to  provide  useful  guidance  in  practical  applications  especially 
if  the  influence  of  the  boundary  conditions  can  be  separately  investigated 
and  if  the  linearized  model  for  the  difference  relations  at  the  interior 
points  are  properly  selected.  Such  fortunate  circumstances  are,  however, 
not  to  be  presumed  in  complicated  equations  systems. 


3.3  Application  to  Navier  Stokes  Equations 

The  Navier-Stokes  equations  system  is  quasi-linear  due  to,  (ij  the 
variable  convective  velocity,  (ii)  the  variable  density  and  energy  and 
hence  the  variable  diffusivity  as  illustrated  in  Equation  (3.3).  It  is 


"l 

ujl 


vAl/M  ~ 075 
Si/jjAf/fAx/-1 » .005 


Fig.  1.  Running-wave  solutions  of  the  non-linear  equation  du/dt  - d2(u2)/dz2.  The 
curves  show  the  exact  solution,  given  by  equation  (8.25)  and  the  dots  show  the  solution  cf  the 
difference  equation  (8.27)  with  6 = 0.4  and  with  At  and  Ax  so  chosen  that  rAf/Ax  ■*  0.075  and 
5M«4A//(Ax)a  m 0.005.  The  numbers  on  the  curves  are  cycle  numbers. 
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further  complicated  by  the  small  diffusivity  or  large  Reynolds  number  and 

f 

the  presence  of  many  such  terms  especially  in  multi-dimensional  flow  prob- 
lems. If  the  standard  procedure  of  local  linearization  is  followed,  the 
resulting  linearized  equations  are  very  long.  The  v. Neumann  stability 
analysis  for  such  equations  will  inevitably  lead  to  unwieldy  algebraic 
expressions  so  that  the  explicit  limit  on  At  at  each  mesh  point  can  only 
be  obtained  at  much  more  labor  than  what  is  required  in  the  situation  of 
Equation  (3.11).  It  is  inqpractical  to  consider  checking  the  stability 
limit  at  many  mesh  points  even  if  infrequently.  It  appears  imperative 
to  look  for  simple  but  meanihgful  model  equations  such  as  Equation  (3.10) 
in  the  previous  example.  This  is  considerably  complicated  by  the  change 
of  the  asymptotic  behavior  of  tho  Navier-Stokes  equations  system  in  differ- 
ent regions  in  the  field  of  computation,  enumerated  in  the  previous  chapter. 
Near  tho  solid  boundaries  or  where  the  viscous  effects  are  important,  the 
region  is  locally  parabolic  or  elliptic.  Far  away  from  the  solid  boun- 
daries, the  direct  viscous  effect  is  negligible  and  the  flow  region  is  pri- 
marily hyperbolic.  It  is  unfortunate  that  a difference  algorithm,  when 
applied  to  practical  differential  equations  of  different  types,  will  lead 
to  difference  equations  with  quite  different  stability  behavior. 

Tables  I and  II  list  a few  common  difference  algorithms  when  applied 
to  the  simple  wave  equation  and  to  the  simple  diffusion  equation  respec- 
tively. An  algorithm  often  yields  a stable  difference  equation  for  the 
diffusion  equation  such  as  the  forward  time-centered  space  algorithm 
(Scheme  1 in  Table  II)  while  it  provides  an  unstable  difference  equation 
(Scheme  3 in  Table  I)  from  the  simple  wave  equation.  Friedrichs  modifi- 
cation, which  renders  the  wave  problem  stable  (Scheme  4 in  Table  1),  on 
the  other  hand,  leads  to  an  instable  diffusion  problem  (Scheme  2 in 
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Table  II).  The  centered  time  and  centered  space  algorithm  is  another 
example  which  is  given  in  these  tables.  There  arc  many  other  examples 
like  these.  Such  schemes  are  therefore  not  useful  for  integrating 
Navier-Stokes  equations. 

There  are  many  other  schemes  which  are  stable  for  both  types  of 
equations  but  At  are  subject  to  different  restrictions  in  different  re- 
gions. Usually  c g—  <1  for  the  wave  equation  and  £ some  fractional 

constant  g for  the  diffusion  equation;  such  as  the  forward  time,  backward 
space.  Scheme  1 and  2 respectively  in  the  two  tables.  The  condition 
c yg  <1  is  known  as  the  Courant-Friedrich-Levy  (CFL)  condition  of  zone 
of  dependence  to  be  satisfied  generally  for  difference  ferns  for  wave 
equations J*^It  states  that  the  zone  of  dependence  of  the  difference  formu- 
lation oust  include  the  zone  of  dependence  of  the  differential  equation. 

When  such  a scheme  is  used  in  integrating  the  Navier-Stokes  equations, 

coiqputational  stability  might  be  expected  if  At  is  locally  chosen  to  be 

[21 

the  more  restrictive  of  the  diffusion  limit  and  of  the  wave  limit: 

At  < ln.flT  • Y ^r]  (3.12) 

where  c is  related  to  the  local  wave  speed,  v is  the  local  kinematic  vis- 
cosity coefficient,  and  y is  some  constant  less  than  unity.  The  precise 
values  of  c and  y azy  be  determined  from  the  v. Neumann  stability  analysis 
of  the  linearized  Navier-Stokes  equations  after  dropping  the  viscous 
terms  or  the  dynamic  terms  respectively.  The  most  restrictive  of  these 
local  limits  over  all  the  mesh  points  may  then  be  taken  as  the  At  for  the 
next  time  increment.  In  actual  computations , it  is  often  necessary  to 
reduce  this  most  restrictive  limit  on  At  further  by  introducing  an  empir- 
ical safety  factor  which  may  have  to  be  rather  small. 
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Some  safety  factor  may  have  to  be  needed  because  of  the  unknown 
effect  of  the  boundary  conditions.  But  actual  computation  often  indi- 
cates the  instability  to  be  initiated  from  the  interior.  Thus,  it 
appears  to  be,  at  least  in  part,  due  to  the  fact  that  the  pure  diffusion 
end/or  the  pure  wave  equation  are  rather  poor  models  for  the  interior 
points  of  the  linearized  Navier-Stokes  equations.  It  is  true  that  in 
the  linearized  form  the  Navier-Stokes  equations  may  be  visualized  as 
the  superposition  of  a wave  and  a diffusion  equation.  The  stability 
limit  is,  however,  not  generally  the  superposed  diffusion  and  the  wave 
limit.  This  is  because  the  determination  of  the  eigen-values  of  a 
linear  equation  is  not  a linear  problem,  which  requires  the  solution  for 
the  roots  of  a polynomial  equation  with  constant  coefficients.  A small 
perturbation  on  a coefficient  of  the  characteristic  polynomial  often  leads 
to  an  ^proportionately  large  change  in  the  largest  eigenvalue  (or  the 
spectral  radius)  depending  on  the  specific  difference  algorithm. 

To  illustrate  the  situation,  consider  the  one-dimensional  Burgers' 
equation  with  constant  c and  v . 


3u 

Tt 


♦ c 


3u  , 
3x 


32u 

nr 


(3.13) 


When  c and  v are  taken  as  the  local  values  at  a mesh  point,  it  serves  as 
a linearized  model  of  the  Navier-Stokes  equation  in  one  space  dimension, 
with  the  essential  characteristics  of  changing  type  of  the  partial  differ 
ential  equation.  If  Equation  (3.13)  is  discretized  with  forward  tine 
backward  space  difference  for  the  convective  term  and  the  centered  space 
difference  for  the  diffusion  term,  the  v. Neumann  stability  limit  is 

At  < (c/L%  ♦ 2v/Ax2V! 


(3.14) 
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which  is  almost  half  of  the  h/perbolic  limit  Ax/c  or  the  diffusion  limit 
Ax2/2v  if  they  are  approximately  equal.  Thus  the  safety  factor  to  be 
applied  to  the  condition  (3.12)  should  be  about  1/2  or  less  for  the 
stable  computation  of  the  interior  points  alone. 

The  situation  is  even  more  critical  for  multidimensional  flow  prob- 
lems. Consider  the  two  different  models  for  2-D  problems 

3u  r3u  * \jr3*u  * 32ui 

3t  + c " V(5*T*  3y2T> 

3u  3u  _ ,.,2* u 32u 

and  3t  C3T  ' v('3*r+  l§p‘)  I 

with  the  convective  term  v in  the  dynamic  terms  represented  by  c 
and  zero  respectively.  The  stability  limits  for  the  two  cases  assuming 
Ax  = Ay  are 


At  £ j (c/Ax  «•  2V/Ax2)-1 

and  At  £ (c/Ax  + 4v/Ax2)-1  ( 

respectively.  This  means  that  the  safety  factor  to  be  applied  to  con- 
dition (3.12)  should  be  M/4  and  1/3  for  2D  flow  problems  and  even 
smaller  for  the  3D  flow  problems.  Thus,  the  stability  condition  (3.12), 
based  on  superposing  the  wave  end  the  diffusion  parts  of  the  Navier- 
Stokes  is  not  very  useful  although  it  is  simple  and  convenient. 

The  local  stability  condition  based  on  the  linearized  Burgers' 
equation  was  found  to  be  quite  satisfactory  for  the  integration  not  only 
of  the  nonlinear  Burgers'  equation  without  a safety  factor,  but  of  the 
full  Navier-Stokcs  equations  when  the  boundary  conditions  are  properly 
treatedP^The  one -dimensional  Burgers'  model  should  be  applied  locally 
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to  the  flow  along  the  streamline  through  a mesh  point.  This  will  yield 
stability  limits  of  the  form  of  Equation  (3.14)  and  (3.16)  in  which  c 
should  be  interpreted  as  the  local  signal  speed  |q|  ♦ |a|.  Here  q is 
the  stream  velocity.  Both  the  local  speed  of  sound  a and  the  kinematic 
viscosity  coefficient  v should  be  evaluated  at  the  local  temperature 
or  energy.  Ax  should  be  evaluated  along  the  streamline  in  some  manner 
and  nay  well  be  taken  as  the  smaller  of  (Ax  , Ay)  for  2-D  problems  for 
example.  The  different  details  how  these  local  quantities  may  be  approx 
imated  by  those  explicitly  calculated  at  each  point,  will  give  large 
varieties  of  expressions  for  evaluating  the  local  stability  limit  on  At 
for  a given  choice  of  the  difference  algorithm  for  discretization.  It 
is  advisable  to  choose  a simple  form  convenient  for  the  explicit  deter- 
mination of  the  limit  on  At  at  each  point,  although  Jess  accurate. 

This  calculation  is  to  be  carried  out  at  many  points  and  at  many  time 
intervals  for  an  estimate  of  the  most  restrictive  limit  (the  smallest 

L 

value  of  ) on  At  for  the  next  time  interval.  It  may  also  be  convenient 
to  check  the  local  linearized  stability  limit  once  every  few  time  steps 
rather  than  at  everytime  step  and  to  adjust  the  magnitude  of  At  adopted 
in  the  confutation  for  the  next  few  steps  accordingly. 

3.4  Treatment  on  the  Boundary 

When  the  appropriate  local  linearized  stability  limit  is  obeyed, 
computational  instability  at  the  interior  points  can  usually  be  avoided  jj 

although  oscillations  of  fairly  large  but  bounded  amplitudes  are  often 
present  in  the  calculated  results.  These  oscillations  originate  from 


the  boundaries,  both  interior  and  exterior  and  do  not  represent  computa- 
tional instability  in  the  sense  of  boundedness  of  the  solution  discussed 
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previously.  Such  bounded  oscillations  are  often  referred  to  as  the 
Nonlinear  Instability  which  is  basically  a different  phenomenon  more 
directly  related  to  the  question  of  accuracy  and  can  doubtfully  be 
clarified  by  the  heuristic  local  linear  stability  treatmen|4Jiscussed 
in  the  previous  section. 

Genuine  unstable  computations  can  result  when  certain  boundary 
treatment  is  applied  to  some  difference  algorithm.  For  such  cases,  the 
local  linearized  analysis  can  often  tell  the  impending  computational 
instability.  Consider  the  integration  of  the  inviscid  gas  dynamic  equa- 
tion (2.7)  with  the  Leap-Frog  scheme,  (Scheme  5 in  Table  I) 


ur1.  iff-1.  -a.  [if iff , l 

) j i s*  Lj*i  j-iJ 


which  is  second  order  accurate  in  both  time  and  space  and  is  always 
stable  for  any  values  of  At/Ax  at  all  the  interior  points.  To  initiate 
tho  integration,  both  U?  and  uj  should  be  available  at  all  j « 0,1,2,...J 
and  boundary  conditions  must  be  provided  at  both  boundaries  j = 0 and 
j * J.  Note  that  both  the  initial  value  uj  and  the  boundary  data  at  j * J 
are  not  specified  by  the  initial  data  of  the  differential  problem  of  the 
propagation  of  a small  wave  in  an  unbounded  flow  field.  These  data  are 
extraneous  and  are  brought  about  by  the  use  of  the  higher  order  accurate 
difference  algorithm  in  which  a first  order  differential  coefficient  is 
replaced  by  a second  order  difference  quotient. 

The  extraneous  initial  data  llj  are  usually  obtained  from  U°  and  the 
temporal  derivatives  through  Taylor  series  about  t * 0.  The  higher  order 
temporal  derivatives  are  obtained  from  the  initial  data  U?  in  the  vicinity 
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of  the  point  based  on  the  differential  relation  and  its  time  derivatives. 

It  is  not  obvious  how  the  extraneous  boundary  data  at  j = J should  be 

defined.  One  of  the  natural  ways  is  to  extrapolate  along  x assuming 
3 

that  -j—  U is  small  is 


>5 


«5-i 


(3.18) 


This  is  not  a bad  physical  approximation.  Computationally  it  leads  to  the 
difference  relation 


«1- 


Si 


-A 


J-l 


At 

Ax 


(3.19) 


for  advancing  the  mesh  value  at  J-l  immediately  preceeding  the  boundary 
point  J.  If  the  v. Neumann  stability  analysis  is  applied  locally  to  this 
difference  equation  with  Aj_^  taken  as  a constant,  and  U taken  as  a scalar 
tnknown,  this  difference  relation  (3.19)  is  locally  always  unstable  with 
the  amplifi cation  factor  |X|  » | IfJ*  L This  is  because  the 

v. Neumann  analysis  leads  to  the  algebraic  relation 

X " X " 'AJ-1  SX  [(*-«>*  **«)  ♦ 1 sin  XAx  J = -2(fr  ♦ if£)  (3.20) 

whero  k is  the  wave  number  under  consideration  and  2ff  and  2f^  are  the 
real  and  the  imaginary  parts  of  the  right  hand  side.  Thus 

A - -(fr  ♦ if*)  ± [l  ♦ (£r  ♦ if^*]172  (3.21) 

For  some  choice  of  k,  f^  will  be  zero  and  the  absolute  value  |k|  will 
be  greater  than  unity  regardless  of  the  magnitudes  of  Aj^  ^ . Actual 
confutation  confirms  the  instability  that  |l/J|  diverges  as  n.  If  Aj_^  ^ 
should  be  taken  as  unity  and  if  the  initial  data  satisfies  » (-l)-*+n 
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for  n * 0 and  1 and  all  j * 0,  1...J-1,  the  solution  of  the  difference 
equation  actually  can  be  shown  to  continue  as 

U?  = (-l)j+n  ♦ (-l)j  F(j  4 n) 


with  F(j  4 n < J)  = o and 

F(n  4 J)  « (-l)n-12n 


Higher  order  accurate  extrapolation  formulas  based  on  zero  higher  order 
derivatives  instead  of  Equation  (3.18)  will  only  change  f and  f^  and 
still  lead  to  computational  instability  in  the  same  manner. 

Careful  examination  of  the  local  stability  analysis  will  suggest 
that  stable  computation  will  result  if  the  extraneous  boundary  value 
l/J  is  obtained  as: 


«5 


1 

7 


fu"*1 


uJ-iJ 


i.e.,  in  the  first  order  extrapolation  formula  (3.18)  is  replaced 
by  the  average  of  its  temporal  neighbors.  Then  the  difference  relation 
on  tho  boundary  is 


i-q 

l+o 


2a 

1+q 


whore  a 


1 . At 
7 AJ-1  Sx 


with  a > 0. 


i *5-1 1 — «*■  ti«5:}i. 

and  the  advanced  values  remains  bounded.  Alternatively  if  the 

v. Neumann  analysis  is  followed,  then 


(3.22) 


(3.23) 


(3.24) 


(3.25) 
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. 1-a  1 2a  -ikAx 

x ‘ Fa  X*  FS  e 

and 

M i Ig  1 fa  * ill 

Thus 

-1  ‘ ■ ps < M < i 

if  a<l 

or 

-1  < lxl  < 1 

if  a>l 

and  computational  stability  can  be  ejected. 

Thus  the  local  linear  stability  analysis  will  help  to  avoid  unfortu- 
nate choices  of  the  unstable  boundary  conditions  and  sometimes  suggest 
appropriate  choices  to  secure  stable  computation.  It  must  be  cautioned 
that  if  the  particular  choice  of  the  boundary  condition  fails  to  represent 
the  physical  situation,  the  computed  stable  solutions  need  not  be  a good 
approximation  to  the  solution  of  the  physical  problem.  It  is  common  that 
oscillations  of  finite  amplitudes  appear  to  be  generated  at  the  various 
boundaries  in  a computed  stable  solution  and  that  such  oscillations  appear 
to  propagate  into  the  interior  of  the  field  of  confutation,  (or  away  from 
the  shock  wave  or  other  interior  boundaries).  They  represent  oscillatory 
error  conf or.snts , superposed  on  the  correct  solution  of  the  physical  prob- 
lem and  are  likely  introduced  by  the  "eirors"  in  the  difference  treatments 
of  such  boundaries.  Indeed  there  are  also  nonos dilatory  errors  caused 
by  the  difference  treatments  on  the  boundary  and  such  errors  may  actually 
be  more  serious  because  of  its  deceptively  smooth  appearance  in  the  re- 
sults of  a stable  confutation.  Such  errors  tend  to  be  overlooked  especially 
in  view  of  the  difficulties  in  securing  a stable  confutation.  An  isfortant 
aspect  of  studying  the  accuracy  of  the  confuted  results  is  to  recognize  if 
the  various  boundary  conditions  are  appropriate  and  to  estimate  the  associ- 


ated errors. 


IV.  IMPLICIT  COMPUTATION  AND  RATE  OF  CONVERGENCE 


Implicit  difference  algorithms  generally  lead  to  stable  difference 
equations  when  applied  to  simple  wave  end  ciaple  diffusion  equations  as  is 
Indicated  In  Tables  I and  II.  The  local  linear  stability  analysis  for 
equation  (3.3)  illustrates  further  that  stability  is  "improved"  when  the 
fraction  6 of  the  spatial  derivative  evaluated  at  the  advanced  tine  level 
(and  hence  implicit)  is  increased  fro®  zero  to  1/2.  The  system  becomes 
unconditionally  stable  when  6 J>  1/2.  The  lepllclt  difference  algorithms 
are  traditionally  used  in  the  solution  of  Laplace  or  Poisson  equations 
without  any  problem  of  computational  stability.  The  implicit  difference 
algorithm  then  appears  to  be  the  most  desirable  from  the  point  of  view 
of  avoiding  computational  stability  especially  for  complicated  problems 
with  mixed  behavior.  It  will  be  demonstrated  that  the  merit  of  the 
implicit  schemes  is  not  really  that  obvious.  There  are  indeed  other 
difficulties  which  may  be  more  serious  than  computational  stability. 

In  implicit  difference  algorithms,  the  difference  relation  at  a given 
mesh  point  contains  the  unknown  advanced  values  of  quantities  at  it  neigh- 
boring nesh  points.  It  is  necessary  to  treat  the  system  of  difference 
relations  at  all  the  interior  nesh  points  simultaneously  to  solve  for  all 
the  unknowns  at  all  the  Interior  mesh  points.  Thus  the  difference  formu- 
lation based  on  a totally  implicit  scheme  will  require  the  inversion  of 
matrices  of  very  large  dimensions.  This  imposes  severe  requirements  on 
the  memory  capacity  and  on  the  arithmetic  speed  of  the  computer.  This 
also  calls  for  skills  in  rendering  efficient  inversion  of  the  sparse  but 
large  matrices  and  inevitably  through  some  form  of  iterative  procedures. 
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The  rate  of  convergence  or  the  number  of  Iterations  required  to  solve  the 
system  of  equations  to  a prescribed  accuracy  Is  of  great  concern.  This  Is 
because  the  computational  effort  required  to  complete  a "sweep"  over  the 
field  of  computation  (l.e.  to  advance  the  values  of  the  functions  at  all 
the  mesh  points  for  one  time  step)  is  generally  much  larger  for  the  im- 
plicit difference  formulation  than  for  the  explicit  difference  formulation. 
It  16  hcped , however,  thet.  In  the  absence  of  a stability  licit  on  At  with 
the  Implicit  difference  formulation,  the  time  steps  may  be  taken  so  cany 
times  larger  than  the  time  step  allowed  by  the  stability  limit  of  the  ex- 
plicit formulation  as  to  more  than  compensate  for  the  much  larger  computa- 
tional effort  per  time  step  for  the  Implicit  formulation.  In  the  following 
sections  this  question  will  be  examined. 

4.1.  Simple  Time  Dependent  Problem. 

The  advantage  of  the  implicit  formulation  Is  best  Illustrated  In  the 
solution  of  the  time  dependent  heat  transfer  problems  In  multispace  dimen- 
sion or  In  the  solution  of  Laplace  equations  for  the  steady  state  problem. 


-§~-  W2u  (4.1) 

?or  such  problems,  the  system  of  simultaneous  difference  equations  to  be 
solved  can  be  conveniently  arranged  to  be 

AD  - f (4.2) 

where  U Is  the  vector  unknown  representing  the  temperature  at  all  tha  N 
interior  mesh  points,  arranged  In  some  appropriate  order,  f is  a known 
vector  of  dimension  H and  A is  an  N x N trldlagonal  matrix  often  diagonally 
dominant.  The  solution  of  the  system  (4.2)  for  the  unknown  vector  U Is 
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Tne  mathematical  foundation  and  the  various  practical  aspects  of  the  numerical  5 
solution  of  gas  dynamic  equations  are  critically  reviewed  with  emphasis  on  obtaining  I 
quantitatively  accurate  solutions  for  application  in  various  engineering  and  sciences  ■ 
Computational  stability  rate  of  convergence  and  accuracy  (or  error  estimate)  are  dis- 
cussed. The  promises  and  problems  of  the  4th  generation  computers  are  outlined  with- 
in this  perspective., 

/C" 

~ Computational  stability  should  not  be  obtained  at  the  sacrifice  of  the  conver- 
gence rate  to  and  the  accuracy  of  the  final  solution.  With  accuracy  in  mind,  the 
explicit  algorithms  are  likely  preferrable  to  the  implicit  ones.  Strict  conserva- 
tion of  the  difference  formulation  is  recommended  and  exemplified  to  avoid  the  accum- 
ulation of  local  trm cation  errors  and  to  facilitate  the  estimate  of  the  errors  in  a 
steady  state  solution.  Illustrative  examples  are  given  including  supersonic  flows 
with  shocks . 
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equivalent  to  the  inversion  of  the  matrix  A,  giving  U as  U • A r. 
Computationally,  a highly  efficient  method  can  be  used  to  solve  the  system 
(4.2)  with  approximately  SN  operation  counts.  (Conventionally,  each  multi- 
plication and  division  counts  as  one  operation  while  addition,  substractlon 
sad  other  data  management  operations  are  neglected).  This  is  to  be  compared 
with  N counts  for  the  solution  of  an  explicit  system.  (The  evaluation  of 
coefficients  is  ignored  here  on  the  presumption  that  the  eaae  amount  of 
computational  work  is  needed  in  both  the  implicit  and  the  explicit  cases) . 
Thus  the  computational  effort  to  advance  the  solution  for  one  time  step 
with  the  implicit  format  is  ebout  5 times  as  much  as  that  with  an  explicit 
format.  As  is  illustrated  in  Table  2,  most  stable  explicit  schemes  will 

possess  the  stability  limit  (easily  verified  by  v.  Neumann  analysis)  of  the 

vAt  'v 

type  a ■ 5—  < (1/2,  1/4);  here  1/4  is  for  two  dimensional  problems  (see 

Ax 

equation  3.16  with  c ■ o.).  With  good  spatial  resolution,  l.e.  a small  Ax, 

1 2 

the  time  step  for  the  explicit  scheme  will  be  limited  to  At  < %"-Ax  which 

is  Indeed  very  small.  Thus,  if  computations  with  the  implicit  formulation 

5 2 

should  be  carried  out  with  a time  ctep  larger  than  — — Ax  or  even  taken 

ee  At  • Ax,  considerable  saving  in  the  computational  effort  results  in  the 

determination  of  the  temperature  field  U at  that  later  time. 

The  benefit  that  results  is,  however,  Illusory  if  the  determination 

of  the  solution  at  cose  cpeciflc  later  time  is  required  to  possets  a specific 

accuracy.  Suppose  that  all  the  variables  are  properly  non-d  linens  ionalized 

-2 

and  that  it  is  required  to  achieve  an  accuracy  of  10  , assumed  to  be  solely 

dependent  on  the  truncation  error  (i.e.  all  the  other  errors  are  suppressed 
In  the  formulation  and  computation).  Then  if  the  explicit  algorithm  1 in 
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table II  is  used,  which  ie  first  order  accurate  in  time  and  second  order 

2 

accurate  in  space  (i.e.  - 0(At,Ax  ),the  field  of  computation  defined 

by  x ■ 0 to  1 and  y ■ 0 to  1 for  a tuo  dimensional  problem  from  time  t • 0 

to  1 should  be  divided  into  at  least  10  equal  parts  in  both  the  x and  the 

y direction,  i.e.  Ax  ■ Ay  » 1/10,  preferably  say  with  Ax  - Ay  “ 1/20  to 

allow  soon  margin  of  safety.  The  stability  limit  will  require  a At  (with 

time  non-diaensionalized  by  the  square  of  characteristic  length  divided 

by  the  diffusivity)  as  spall  as  — y Ax^  - 10  ^ if  Ax  ■ 1/20  or  1/4  x 10  ^ 

if  Ax  ■ 10  The  relative  magnitudes  of  At  and  Ax  are  such  that  the  local 

2 

truncation  error  efc  • 0(At,Ax  ) end  hopefully  the  computed  results  will 
remain  consistent  with  the  accuracy  requirement.  (In  this  case,  the 
accumulation  of  the  local  truncation  error  will  remain  of  the  same  order.) 

If  now,  the  scheme  (1)  in  Table H is  modified  so  that  the  spatial 
derivative  is  replaced  by  the  implicit  difference 

•<  > 

2 

for  both  x and  y direction  with  the  same  local  truncation  error  efc  ■ 0(At,Ax  ). 

This  schema  (Loasonen)  is  unconditionally  stable,  i.e.  At  can  be  taken 

arbitrarily  large  compared  vithAx  without  suffering  computational  instability. 

2 

However,  with  At  much  larger  than  Ax  , the  local  truncation  error  is  of 

0(At)  » O(Ax^).  Thu 8 with  Ax  * 1/20  as  was  in  the  explicit  case  and  if 

At  Is  taken  as  Ax/5,  which  is  16  times  larger  than  the  stability  limit  of 

the  explicit  scheme,  the  computational  effort  will  be  only  1/3  of  that 

with  the  explicit  scheme.  The  solution  so  obtained  is;  however,  less 

accurate  with  e£  * G(At)  and  At  * -y—  • 10  , marginally  acceptable  to 

_2 

the  required  accuracy  10  , allowing  no  room  for  the  accumulation  of  the 

local  truncation  errors.  Formally,  this  solution  from  the  Implicit  scheme 
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should  be  compared  with  the  solution  from  the  explicit  scheme  with  Ax  ■ 10  * 
2 2-2 

with  e^  ■ 0(Ax  ) and  Ax  ■ 10  . The  computational  effort  of  this  explicit 

scheme  is  actually  only  80Z  of  the  Implicit  scheme  with  the  came  local 

truncation  error.  Alternatively  if  the  implicit  scheme  is  to  produce  a 

result  with  accuracy  comparable  to  the  explicit  solution  computed  with 

1 2 

Ax  - 1/20  and  At  « 1/4  x(  ) . The  time  step  At  for  the  implicit 

2 1 2 

calculation  should  be  taken  at  most  as  Ax  - ^q7J~  that  et  - 0(At  ■ Ax  ). 
Then  the  computational  effort  for  the  explicit  format  will  again  be  80Z 
of  that  of  tho  implicit  format  of  comparable  accuracy. 

The  effectiveness  of  the  implicit  algorithm  is  largely  nullified  by 
the  first  order  temporal  accuracy  of  the  difference  scheme  in  the  above 
example.  It  may  be  that  implicit  schemes  with  second  order  temporal 
accuracy  will  be  more  effective  in  reducing  the  overall  computational 
effort • but  such  higher  order  schemes  will  be  cumbersome.  Prom  this 
point  of  view  alone,  the  implicit  schemes  would  appear  to  be  certainly  ad- 
vantageous in  the  solution  of  steady  state  problems  via  asymptotic  temporal 
approach  since  the  temporal  accuracy  is  of  little  concern  But,  as  will  be 
discussed  in  the  next  section,  it  is  not  certain  if  the  large  temporal  steps 
is  conducive  to  rapid  convergence  to  the  steady  stace.  It  should  be  noted 
in  the  above  example,  that  the  solution  of  the  implicit  formulation  calls 
only  for  the  inversion  of  a trldlagonal  matrix  which  can  be  implemented 
most  efficiently  in  o,  5N  operations.  For  fluid  dynamics  problems,  the 
matrices  resulting  from  an  implicit  formulation  will  be  far  more  complex 
and  tha.  solution  of  such  matrix  equations  will  be  far  more  time  consuming. 

It  appears  prudent  not  to  expect  significant  aavings  in  the  computational 
effort  by  the  use  of  implicit  difference  algorithms  without  some  detailed 
investigation. 
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4.2  Iterative  Solution  of  Steady  State  and  Asymptotic  Temporal  Approach 

Most  of  the  fluid  flow  problems  of  practical  Interest  are  at  steady 
state  or  quasi-steady  state  In  which  the  temporal  variations  of  the  flow 
variables  are  negligible.  Discretlzaticn  of  such  steady  state  equations 
will  generally  lead  to  implicit  difference  relations  in  terms  of  the  steady 
state  values  of  various  phyoical  quantities  at  all  the  interior  points  and 
the  boundary  values.  Except  for  the  solution  of  potential  flow  problems 
of  Incompressible  fluids,  the  differential  equations  will  be  non-linear 
and  considerably  more  complicated  than  the  Laplace  equation.  The  re- 
sulting implicit  difference  relations  will  give  rise  to  a rather  sparse 
matrix  A,  when  written  in  the  format  of  equation  (4.2).  The  sparse  matrix 
A will  not  be  trldlagoaal  or  block-tridiagonel,  or  other  spcial  forma  convenient 
for  the  solution  of  the  system  of  equations.  In  fact  the  nonlinear  terms 
will  first  have  to  be  quasi-linearized  so  that  the  coefficients  in  matrix 
A can  be  evaluated  with  some  assumed  approximate  values.  The  system  of 
linear  equations  will  then  be  solved  iteratively  until  the  solution  from 
(4.2)  agrees  with  the  assumed  solution  under  certain  convergence  criteria. 

Let  superscript  n indicate  quantities  evaluated  with  the  n*'*1  iterate 
of  U and  use  the  system  of  difference  relations  (4.2)  to  calculate  the 
(n+lj^ iterate . Then  equation  (4.2)  becomes 


*V+l  - £°  . 0 


(A. 3a) 


which  is  Indeed  the  same  as 


A0^1-  U°)  - fn  - AV1 


(4.3b) 


Equation  (4.3b)  can  now  be  considered  as  obtained  from  a tine  dependent 
equation  In  which  terns  with  spatial  derivatives  are  the  sane  as  those  in 
the  otcady  state  equation  (4.2)  but  with  an  added  temporal  term 


- SO  - 


lim  ,.,ii  Un+1-Un  v 3u 

At  - O AtA  ~Tt — ■ AtA(u)  TT 


3 u 


with  forward  temporal  difference  quotient  replacing  g-£—  . 

The  iterative  solution  of  a steady  state  problem  based  on  en  implicit 
algorithm  then  is  not  substantially  different  from  the  solution  of  a time 
dependent  problem,  albeit  the  artificial  temporal  term  may  net  correspond 
to  the  temporal  terms  in  the  time  dependent  form  of  the  Ncvier-Stokes 
equations.  The  physical  meaning  of  the  individual  fictitious  temporal 
terms  can  be  easily  identified  when  the  matrix  operator  A is  written  in 
expanded  form.  Thus  the  iterative  index  n can  be  identified  with  the 
temporal  index  n in  the  time  dependent  formulation  although  the  equivalent 
time  dependent  physical  problem  may  contain  artificial  sources  of  mass, 
■omentum  and  energy.  These  artificial  sources  are  small  but  distributed 
over  the  entire  field  of  computation, In  the  interior  as  well  as  on  the 
boundary,  and  vanish  in  the  steady  state  limit. 

In  the  numerical  solution  of  the  Navler-Stokes  equations  in  multi- 
space  dimensions,  there  will  be  a few  thousand  mesh  points  and  4 or  5 
unknown  quantities  at  each  mesh  point.  The  dimension  N of  the  vector 

4 

U is  commonly  0(10  ) or  larger.  To  oolveequ&tion  (4.3a)  for  the  successive 
approximations  to  the  solution  of  the  nonlinear  equation  (4.2)  at  each 
time  or  iterative  step  by  the  standard  Gaussian  elimination  process, 
requiring  ± N operations  per  step,  is  out  of  the  question.  It  is 
imperative  to  develop  highly  efficient  iterative  methods.  Thus  in  equations 


- 51 


(4.3),  the  matrix  operator  An  is  split  into  two  parts  with  Bn  operating  on 
U°+1  and  (An-Bn)  operating  on  U11 . Thus 


BnUn+1  + <An-B°)  U“  - fn 


(4.4a) 


Bn(Un+1  - u°)  - fn  - aV1  (4.4b) 

where  Bn  may  be  some  conveniently  invertible  matrix  so  that  Un+'^  can  be 
conveniently  solved.  This  un+^  will  replace  Un  in  the  next  Iteration  until 
finally  Un+^i  0°  according  to  some  steady  state  criterion.  In  this  manner, 
the  Iterative  solution  of  the  qucsA-linearized  equation  (4.3)  has  incor- 
porated the  iterations  that  were  called  for  by  the  quasi-linearization  of 
the  nonlinear  equation . 

If  B is  chosen  ao  the  identity  matrix  I,  equation  (4.4b)  becomes 
identical  with  the  difference  equation  obtained  from  the  explicit  scheme, 
using  the  forward  time  difference  algorithm  and  the  spatial  difference 
algorithm  of  the  implicit  equation  (4.3a).  Thus,  like  the  solution  of 
time  dependent  equations  with  explicit  schemes,  the  Iterative  solution 
of  equation  (4.4b)  with  B ■ 1 corresponds  to  tracing  the  physical 
development  of  the  time  dependent  flow  field  from  an  Initial  state  toward 
the  oteady  state.  The  local  accumulations  of  mass,  momentum  and  energy 
in  the  cell  around  each  mesh  point  are  precisely  how  they  would  be 
expressed  in  the  explicit  scheme  for  the  time  dependent  flows. 

The  matrix  B may  be  chosen  to  be  diagonal  with  its  diagonal  elements 
equal  to  the  diagonal  elements  of  A,  l.e.  b^  “ and  b^  * 0 for  i f J. 
Such  an  iterative  process  is  known  ao  Jacobi  iteration.  Since  b^  - a^ 
are  not  identically  unity,  the  temporal  terms  may  be  larger  (or  smaller) 
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than  the  accumulation  term  in  the  physical,  time  dependent  flow.  The 
excess  (or  deficiency)  of  that  particular  quantity  U may  be  attributed  to 
the  presence  of  a source  (or  sink)  of  that  quantity  at  the  mesh  point  under 
consideration.  These  artificial  sources  (sinks)  will  tend  to  zero  when 
the  asymptotic  steady  state  is  approached. 

If  the  matrix  B is  chosen  to  be  tha  main  tridiagonal  elements  of  A,  l.e. 
b^  ■ for  |i-j|  _<  1 and  b^  ■ 0 for  |i-j|  > 1;  then  the  artificial 

temporal  terms  will  contain  spatial  derivatives.  They  represent  the 
doublets  and  quadruples  of  the  source-sink  pairs  around  the  mesh  point. 

The  situation  is  quite  complicated  algebraic ially  and  physically,  but  it 
is  vary  natural  physically  how  a steady  state  may  be  reached  via  such 
time  dependent  states  provided  that  all  these  sources  and  doublets  etc. 
properly  vanish  in  the  steady  state  limit.  In  practice  the  choice  of  B is 
dictated  by  the  desire  to  reduce  tha  computational  effort  in  obtaining 
the  steady  solution.  Irrespective  of  its  physical  correspondence  to  some 
temporal  flow  field.  The  purpoce  here  is  to  show  that  the  asymptotic 
temporal  approach  end  the  iterative  solution  of  the  Implicit  formulation 
to  obtain  steady  state  results  are  fundamentally  similar.  The  Iterative 
method  does  take  much  larger  computational  effort  per  iteration  or  per  time 
step.  But  it  permits  the  use  of  a much  wider  variety  of  temporal  arti- 
fices to  produce  a very  rapid  convergence  to  the  steady  state,  possibly 
with  less  overall  computational  effort.  It  is  possible,  of  course,  that 
for  some  choices  of  B,  there  may  not  be  any  steady  state  solutions  or 
there  may  be  steady  state  solutions  different  from  what  is  desired  as 
the  corresponding  physical  situation  may  suggest. 


i 
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4.3  Iterative  Methods 

One  of  the  most  popular  choices  of  D is  the  lover  triangular  matrix 
part  of  A.  i.e.  bjk  ■ 0 if  k > J and  - a^  if  k j.  ThJ.3  is  the 
G&uss-Seldel  iteration  or  successive  relaxation  procedure.  The  (n+l)t*1 
iterate  is  given  as 


-JL.  ( f“  _ K1  ijO+1  _ ? n n"  x 

ajj  1 tSl  ik  k U Jk  k 


(4.5a) 


from  which  the  successive  scalar  components  of  TJ11*1  can  be  explicitly 
calculated  in  the  order  of  increasing  j in  which  the  latest  available  mesh 
values  are  used  throughout.  This  semi-explicit  solution  of  U11^1  can  be 
given  in  matrix  fora  aa: 


o"+1  - if  + (BV1  <f°  - aY) 


(4.5b) 


Here  (f n-An0n ) is  the  residue  and  (Bn)  ^ is  the  Inverse  of  the  matrix  Bn. 

If  the  vector  calculated  from  equations  (4.5)  is  taken  as  a provisional 
solution  and  if  the  new  Iterate  U°+^  is  evaluated  as  ocme  weighted  average 
of  l end  this  provisional  value  with  weights  (1-0 ) and  0 respectively. 


then. 


' “1 

on+1  - (1-6)  0°  + 610°  + ( Bn)  (f“  - aV1)] 


(4.6a) 


which  is  the  same  as 


0°+1  - U°  + 6(B°f  (fn-AnU°) 


(4.6b) 


or  - u“  + (Bn/6)~1(fn-Anll0) 

The  0 is  often  called  the  acceleration  or  relaxation  parameter.  Equation  (4.6b) 
suggests  that  6 may  be  interpreted  alternatively  as  a multiplier  of  tho 


residue  In  equation  (4.5)  or  of  the  operator  Bn  of  the  artificial  temporal 
sources  in  equation  (4.4b).  The  purpose  Is  to  effect  a faster  convergence 
of  the  Iterative  sequence  by  choosing  an  appropriate  value  of  3.  This 
process  Is  referred  to  as  successive  over  (or  under)  relaxation  when  3 > 1 
(or  3 < 1).  For  the  integration  of  Laplace  equation  in  a rectangular 
domain,  the  optimal  relaxation  parameter  3*  for  the  fastest  convergence 
can  be  evaluated  from  the  mash  spacing  and  Is  usually  around  1.8  - 1.5. 

For  more  complicated  situations,  the  choice  will  have  to  be  empirical 
and  the  optimal  choice  need  not  even  be  an  over-relaxation.  Unfortunate 
choices  of  3 could  lead  to  diverging  sequences,  even  for  Laplace 
equation  (l.e.  beyond  2 >^3^0). 

2 

Each  Gauss- Seidel  Iterative  step  rsqulsss  H operations.  This  Is  to 

3 

be  compared  with  the  counts  of  N /3  for  Gauss  elimination  solution  for  the 
quasi-llnear  steady  state.  The  iterative  solution  would  be  advantageous 
If  It  converges  within  N/3  iterations  since  the  nonlinear  Iterations 
for  the  solution  of  equation  (4.3)  would  then  be  avoided.  Now  with 
N - 0(10  ),  it  Is  hoped  that  by  proper  choice  of  the  relaxation  parameter 
3,  much  fewer  iterations  than  N/3  may  be  needed  to  reach  a steady  state. 

In  principle,  if  the  steady  state  is  defined  by  ||oB+*'-  l^ll/llu11  ||  < 10”°, 
the  number  of  iterative  steps  required  to  converge  csn  be  estimated  by 
m/R  where  R Is  the  rate  of  convergence  with  R - logj^(  — ) and  p is  the 
geometric  mean  of  the  spectral  radii  of  the  matrices  (I?1)”1  An  at  successive 
Iterative  steps  n.  Such  an  estimate  of  R is  not  possible  In  practice  be- 
cause of  the  complexities  of  the  matrix  A and  its  dependence  on  the 
solution  U°. 
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n For  the  Integration  of  some  form  of  hydrodynamic  equations.  It  la  not 
uncommon  that  hundreds  of  such  Iterations  were  needed.  This  is  partly  due 
to  the  nonlinearity  of  the  equations  system  and  partly  due  to  neei-optimal 
choices  of  the  relaxation  parameter.  It  la  also  often  that  such  iterations 
fall  to  converge  despite  vide  range  of  choices  of  the  acceleration  parameter* 
Nov,  if  the  physical  state  of  the  flow  is  eteady  or  quad-steady,  the  asym- 
ptotic temporal  approach  of  using  the  correct  time  dependent  equations 
(B»I)  may  be  expected  to  converge  on  purely  Intuitive  grounds,  provided 
that  the  difference  system  Is  stable  and  ccnslotent  with  the  time-dependent 
Nnvier-Stokes  equations.  But  when  the  Implicit  Iterative  method  Is  used. 

Its  convergence  to  the  steady  state  cannot  be  presumed  on  physical  grounds. 
The  artificial  sources  of  mass,  momentum  end  energy  are  Introduced  purely 
algebraically.  The  particular  temporal  variations  of  these  sources  need 
not  provide  any  steady  state,  although  without  such  external  artificial 
sources,  nature  has  demonstrated  that  a steady  state  will  eventually  be 
reached.  It  might  even  be  legitimate  to  question  If  the  steady  state 
so  reached  should  be  the  same  as  the  one  reached  under  zero  external  sources 
since  the  time  Integrals  of  the  artificial  sources  may  have  altered 
appreciably  the  Integrals  of  motion  of  the  system.  It  la  regrettable  that 
no  useful  soever  con  be  derived  physically. 

Mathematically  speaking,  the  matrix  B in  equation  (4.4b)  can  be  quite 
arbitrary  end  chosen  in  a great  many  different  ways  and  even  be  chosen 
differently  for  different  steps.  Convergence  to  steady  state  Is  assured 
provided  that 
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and  [<Bn)”1An]l<Bn“1r1  A15'1)]  . . . [ (iP  )’1(A°)]  -0  (4.7b) 

n -1  n “1  n 

They  con  be  secured  If  the  spectral  radii  of  all  (B  ) and  (B  ) A are 
less  than  unity.  If  the  fora  of  B chosen  should  be  the  sane  for  all 
iterations,  (4.7b)  is  not  really  much  different  from  the  local  linearized 
stability  criterion  of  v.  Neumann  with  the  amplification  matrix  G re- 
placing B_1A.  (See  Chapter  III,  Section  3.1  end  3.2).  The  essential 
difference  lies  then  in  the  freedom  of  choices  of  the  form  of  the  matrices 
Bn,  at  different  Iterative  steps  n.  It  is  not  clear  if  the  condition 
(4.7a)  would  Imply  the  physical  requirement  of  the  conservation  of  the 
integrals  of  motion.  It  is  also  not  practical  to  find  the  spectral 
radii  or  bounds  of  the  eigen  values  of  these  complicated  matrices.  There 
le  no  counterpart  of  the  local,  linear  stability  analysis  to  provide  oome 
idea  of  the  rate  of  convergence  of  a complicated  problem.  There  is  only 
the  practical  solution  of  trying  it  out  on  the  computer. 

In  practice,  the  possible  choices  of  the  form  of  B ie  severely  limited 
by  its  being  easily  invertible  to  facilitate  the  computation.  It  is  diffi- 
cult to  find  one  that  may  show  significant  improvement  over  the  optimal 
overrelaxation  process  If  the  inferences  from  the  study  of  the  solution 
of  Laplace  equation  is  any  guide,  further  substantial  reduction  of  operation 
counts  per  Iteration  step  is  derived  from  cyclic  processes  built  upon  the 
Gaus8-Scidel  Iterative  procedure.  If  the  field  of  computation  constitutes 
p columns  of  q elements  in  each  row  with  p q ■ N,  the  matrix  B may  be 
chosen  block-lower-triangular  so  that  each  of  the  q blocks  consists  only 
the  p (or  q)  elements  in  each  column  (row).  Then  Bn  Bn  . ,B° 

can  be  assigned  as  the  lover  triangular  matrix  in  successive  blocks  with 
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zero  elements  everywhere  else.  This  Is  the  line  Gauss-Seldel  process  oper- 
ating on  successive  columns,  (or  rows).  Such  a line  process  can  be 
accelerated  by  employing  some  proper  acceleration  parameter. 

The  line  processes  along  columns  and  rows  (diagonals  or  other  con- 
venient directions)  may  be  employed  in  succession  such  as  the  sequence  of 
operators  (B°  ^ . . . Bn  <*)(Bn~<*+*.  . . Bn)  and  its  cyclic  repetition.  A 

set  of  acceleration  parameters  may  be  employed  with  the  cyclic  column-row 

sequence.  This  is  known  as  the  alternating  direction  methods  or  the 

[6] 

method  of  Peaceman  and  Rachford  who  first  demonstrated  the  power  of  such 

cyclic  iterative  methods  for  the  solution  cf  Laplace  equations.  Such  line 

methods  derive  the  benefit  of  less  computational  work  from  the  basic  fact 

that  the  operation  counts  of  Gauss-Seldel  process  is  proportional  to  the 

"square"  of  the  vector  length  of  the  unknown.  Thus  the  operational  counts 

1/2 

of  a complete  cycle  is,  with  p • q-  N for  example, 

P*<j2  + q.p2  - (p+q)pq  ^ 2N3^2  (6.8) 

compared  with  the  N2  for  the  point  Gauss-Seldel  process.  This  means  a 

decrease  of  the  operation  counts  per  sweep  by  the  f actor  2 //n"”,  significant 

3 

to  the  order  of  magnitude  with  N ■ 0(10  ).  The  extension  of  such  cyclic 

1/3 

process  to  problems  in  three  space  distensions  with  fA'q'Vz'vtf  is  obvious. 

4/3 

The  total  operational  counts  per  cycle  is  ^ 3N  and  the  factor  of 

2/3 

operational  counts  reduction  will  be  3/N  . The  preference  of  alternating 

direction  iteration  (ADI)  or  any  such  cyclic  line  iteration  process  over 
the  point  Gauss-Seldel  relaxation  process  is  clear.  The  success  of  reducing 
the  overall  computational  effort  in  the  solution  of  steady  flow  problems 
with  such  schemes  requires,  in  addition,  the  appropriate  choice  of  the 


58  - 


acceleration  parameters  suiteble  for  the  type  of  problems  with  the  class  of 
prescribed  boundary  data.  This  is  where  the  uncertainty  resides. 

For  the  solution  of  Laplace  equations,  the  optimal  acceleration  para- 
meters and  the  maximum  rates  of  convergence  of  these  processes  con  be 
explicitly  determined.  The  ADI  process  is  certainly  the  most  efficient. 

The  sane  is  likely  to  be  true  for  the  integration  of  purely  elliptical 
equations  especially  those  with  the  Laplacien  operator  as  the  leading  terms. 
For  more  complicated  equations,  including  the  equations  of  hydrodynamics, 
much  depends  on  the  ability  of  selecting  the  appropriate  acceleration 
parameters  for  the  problem  at  hand  aud  on  how  the  boundary  conditions  are 
implemented.  Fox  hyperbolic  problems  with  discontinuous  solutions  as  interior 
boundaries,  success  irith  the  implicit  methods  is  yet  to  be  demonstrated. 


3 
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4.4  Fractional  Time  end  other  Alternating  Direction  Methods 

An  alternating  direction  iterative  method  has  been  developed  extensively 

[51 

in  the  Soviet  Union  by  Yanenko,  Marchuk,  etc.,  known  as  the  time  splitting 
or  fractional  time  step  methods.  The  key  idea  is  to  split  the  operator 
as  a sum  of  implicit  difference  operators,  each  of  which  will  lead  to  an 
easily  invertible  tridiagonal  matrix.  The  successive  split  operations  in 
a complete  cycle  serves  as  a "weak"  approximation  of  the  original  operator. 
They  prefer  unconditional  computational  stability  and  formal  second  order 
accuracy  of  the  Crank-Nicholaon  algorithm.  This  was  illustrated  in  Equa- 
tion (3.6)  when  9 takes  the  value  y . Second  order  accuracy  is  needed 
since  a first  ordGr  accurate  scheme  can  hardly  meet  the  accuracy  requirement 
of  practical  problems  with  the  currently  available  computing  machine. 

Its  development  and  its  relative  merits  pertinent  to  the  gas  dynamic 
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applications  are  presented  below. 
Consider  first  the  equation 


(4.9) 


where  L is  a linear  spatial  differential  operator  explicitly  independent 
of  time  t.  Discretized  with  the  Crank-Nicholson  algorithm  which  is 
second  order  accurate  in  both  ties  and  space, 


.u+1  .n  An+1 . 

V f 0 

Let  1 be  the  identity  operator,  we  have 


(4.10) 


(I  + -f5-  L)  $°+1  - (I  - ~ L)  «n 


or 


.n+1  . At  . .-I  » T A.  . . .n  _ >n 

♦ - (I  + — j L>  <r  “ I-  L>  <♦>  " C ♦ 


9 ^ At 

For  the  simple  heat  diffusion  equation  L n,  -a  — , the  matrix  (1  + —5 — L ) 

X r.  Jm  4 X 

9X 

is  trldlagcnal  and  the  spoctrsl  radius  of  the  matrix  C is  obtained  as 


P(C) 


1-S  ^ 
1+S  <|> 


with  ip 


) 


when  0 x • JAx  _<  (J+l)Ax 
Thus 


1-S» 

1+S1{» 


V 


- Vl+S4<  * 


I mi 


which  establishes  the  boundedness  and  unconditional  stability.  Indeed 
with  At/Ax  taken  as  constant,  the  computational  error  is  bounded 
||e||  ^ I leQ  II  + 0(Ax2).  When  the  C.F.L.  condition  At/Ax  <_  1 for  the 
wave  equations  is  satisfied,  this  scheme  was  expected  to  work  for  both  the 
diffusion  and  the  wave  equations,  and  was  hoped  to  work  for  the  Navier-Stokes 


■M 
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type  equations  at  least  for  the  1-D  case,  (which  Is  not  necessarily  true 
as  is  shown  in  section  3.2  for  other  difference  algorithms). 

Consider  now  the  heat  diffusion  problem  in  three  space  dimensions 


0. 


New  L - + Ly  + or  + l^  + L^.  While  the  operations  (I  + L^) 

can  be  easily  inverted,  the  combined  matrix  [I  + (L^  + + L^)] 

is  no  longer  tridiagonal,  and  although  highly  sparse,  cannot  be  simply 
inverted.  So  the  equation  is  to  be  integrated  in  three  successive  steps 
for  the  time  Interval  tQ  _<  t tn+3  and  is  formally  designated  as  the 
fractional  tine  steps  t^^.  tn+2/3  end  tn+3/3  - t^.  For  each  step 
at  tQ  + q^3  the  Crank  Nicholson  algorithm, 

^n+a/3  „ (1  + L^-l  (I  £|_  L^) 
is  used,  thus: 


,n+l 


3 

- n 

a-1 


(1  + 


At 


L ) 
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-1 


(I 


At 
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K > * 


- {I  - AtL  + [ L2  + Z l (LL.  - LflL  )+  . . .]  +0(At3)}<J>n 

2 a-1  0-atl  a 3 3 a 

and  Is  -(I  + L)_1(I  - L)  + 0(At3)}  *n 

if  L Lais  ccemutablc. 
a p 

Thus  the  split  difference  scheme  will  be  second  order  accurate  if  the 
split  operators  are  commutable.  Otherwise,  it  is  only  first  order  accurate. 
Whan  such  commutativity  of  split  operators  for  different  dimensions  (x,y  & s) 
is  not  true,  the  split  scheme  of  only  first  order  accuracy  can  give  second 
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order  accurate  results  In  two  cycles  If  the  cycle  Is  repeated  in  the  opposite 
direction.  For  the  two  consecutive  cycles , i.a. 

*n+1  - n (i  +^-  V1  « - -r  La  > *n 
0-1 


.H+2  -pi  /T  . At  T S'ifv  At  T % aP*1 

♦ - ? ,(1  + — V (I  - “ lo  } * 


the  two  non-commutative  terms  would  cancel  and  the  second  order  accuracy 
is  resumed  for  non-coiasutatlve  operators  L^' s . This  statement  will  ba 
true  even  if  L^'s  involve  differential  operators  with  varying  coefficients 
or  even  depend  on  $ for  ga3  dynamic  quasi-linear  equations  so  long  as 
such  coefficients  are  smooth  and  properly  treated. 

For  unconditional  stability,  it  is  required  tbet  the  operator  L 
and  all  the  split  operators  L^,!^  & are  semi-positive  definite,  that  is, 
the  inner  product  of  (L$,$)  0 for  any  arbitrary  function  and  defined 

over  the  entire  field  of  computation.  This  condition  is  crucial  in 
securing  unconditional  computational  stability.  Take  the  norm  of 
and  define  the  norma  of  an  operator  as  the  natural  norm  deduced  from  any 
vector  norm  then 

l (1  + -jr  L)"1  (I-  nb  L>  (I  + %-  - — | L)  ] 

ll^ll2- * II* 

( «n,oa> 


Define 


(i  + l)"1  - r1 


(I  +-^|  L)C°I 


■2  - a2ii*”h2. 
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Here 

1 | (I  - ^§“  LKn|  |2-  H ^ L)£  Q,  (I  - L)  5n] 

- |Un||2  - At  [L(On.5n]  +-|-||L(Cn)  ||2 

| |(I  + ^§-L)  Cn||2  - [(I  + -~LKnt  a +^L)  en] 

- | |£n|  |2  + At[L(5r‘),  5n]  + ^ IlLCC0)!!2. 

2 

This  A corresponds  to  the  pqu^re  cf  the  spectral  radius  p(C)  for  the  simple 

1-D  heat  diffusion  problem  with  the  Crcnk-Nicholson  algorithm.  Since  both 
2 2 

IUII  «d  ||L«)|r  are  positive  and  At  > 0$  the  amplification  factor 
2 

A will  fca  ^ or  < than  unity  depending  on  whether  (L  (?) , £)  >_  or  < 0 

The  successive  application  of  the  split  operator  at  each  ctc.p  leads 
to 

Il♦n+1||2-A1VA^2II*“!I:! 

2 2 2 

The  conditions  for  A ^ ,Aj  & A , to  be  less  than  or  equal  to  unity  are 
equivalent  to  the  seal-positive  definitions  of  the  split  operators  L^, 

Lj  & Lj,  l.e.  (La$,<£)  _>  0 for  a * 1,2,3.  It  is  true  then  that  unconditional 
stability  results  if  all  operators  LQ  of  L are  semi -positive  definite, 
l.e.  ^ 0.  Thio  sexi-poeitive  definiteness  is  a sufficient  condition 

for  stability  for  tho  entire  cycle  but  not  a necessary  cne. 

How  if  this  restriction  of  serai-positive  definiteness  is  to  be  assured 
this  will  practically  limit  its  applicability  to  the  simple  diffusion 
equation  or  the  Laplace  equation  in  rectangular  domain  and  for  the  Dlrichclct 
problems.  It  is  only  a matter  of  intuition  that  the  method  should  be 
applicable  to  a wider  class  of  circumstances  than  those  for  which  proofs 
could  be  given.  The  situation  is  really  no  better  off  than  vhat  was  en- 
countered for  the  explicit  schemes  on  the  question  of  stability.  There 
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io  indeed  not  even  a necessary  criterion  of  computational  stability  of 
this  split  method  comparable  to  the  von  ITeumar.n  stability  criterion  for 
explicit  schemes. 

In  the  theoretical  treatment  of  the  gas  dynamic  problems,  (Marchult, 
Ysosnko,  etc.)  the  seml-poeitivc  definite  condition  ia  satisfied  by  imposing 
the  special  "periodic”  boundary  conditions  to  the  problem;  in  which  case 
it  is  clear  that 

(!*♦»♦)  - 0 
2 2 2 2 

eo  that  A -A^  “^2  "^3  “ the  »om  is  preserved,  i.e. 

Il*n+1ll  - nvil-  • • • ■ ll*°ll 

This  appears  to  be  an  excellent  feature  for  initial  value  problems.  But 
it  also  implies  that  whatever  error  in  the  initial  data  (if  it  is  a guess) 
will  cot  decrease  in  the  mean  square  norm,  for  example.  The  splitting 
scherae  should  not  be  used  to  obtain  steady  state  aolutlocs  with  the  periodic 
boundary  condition  because  the  results  will  never  be  better  than  the 
initial  guess  within  the  integral  norm. 

For  treating  practical  problems  the  physical  significance  of  this 
stability  requirement  of  (L<t>,4)  >,  0 need  ba  more  carefully  examined. 
Poat-multiply  with  <J>  the  equation 

+ 14  ■ 0 

and  sum  (or  integrate)  over  the  entire  field  of  computation  to  give 
■Jt  11*11+  ■ 0 

X • ® e 


JZ  I Ml2  - - 04  A) 


If  L is  semi-positive  definite,  i.e.  (I/J>,<j>)  ^0,  then  ||<J>||  j<  0. 

This  of  course  implies  the  boundedness  of  the  solution  <p  at  all  tires  and 
a decreasing  sequence  of  j|  <£  ||2.  Now  the  gas  dynamic  equations  in  primary 
physical  variables  are  such  conservation  lays  for  mass,  momentum,  and 
energy. L$  is  the  net  fluxes  of  these  conserved  quantities  out  of  unit 
physical  volume.  If  <{>  and  hence  L<$  are  periodic  over  the  parallelepiped 
in  the  physical  space  to  secure  (L$,4>)  • 0,  then  the  outflux  and  the 
Influx  across  the  boundary  of  computation  exactly  balances.  Thus  with 
4>  identified  as  mass*  momentum  and  energy,  this  condition  excludes  the 
loss  of  mass,  momentum,  and  energy  throughout  the  entire  field  of  com- 
putation. This  means  that  tha  computed  results  should  not  be  expected 
to  show  body  forces  acting  on  some  Immersed  body  or  heat  transfer  to 
and  from  the  body.  If  there  should  be  any  in  the  computed  results  with 
periodic  boundary  conditions,  whatever  lift,  drag  and  heat  transfer  as 
may  be  present  in  the  computed  reaults  originate  from  some  computational 
artifices  and  are  physically  meaningless. 

If  now  the  flow  field  Is  computed  with  periodic  (asymmetric)  boundary 

2 2 

conditions  in  the  transverse  plane,  then  A2  - A^  ■ 1.  A deficit  of 

2 

the  out-flux  14  sad  a positive  $ will  render  (L^,<j>)  j<  0 and  hence  A^  > 1 

2 2 

when  there  is  a body  dmg  (or  energy  sink  to  the  body) . Hence  A - A^  > 1 

ar.d  the  computation  will  be  "unstable".  To  secure  computational  stability 

under  tha  circumstances,  it  is  necessary  to  modify  the  boundary  conditions 

2 2 

in  the  transverse  plane  so  that  A,,  & A.  ^ are  sufficiently  smaller  than  unity 

2 2 2 

to  render  A^  A2  A ^ <1.  Guidance  is  badly  needed  here  for  handling 

these  boundary  conditions  properly  to  secure  computational  stability  with 
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the  split  schemes.  And  even  if  computational  stability  is  achieved  by 
some  means,  there  is  little  idea  how  the  results  of  such  important  quantities 
as  body  drag,  lift,  and  heat  transfer,  as  calculated,  will  compare  with  the 
physical  situation. 

The  fractional  time  step  method  outlined  above  was  not  meant  to  be 
applied  to  the  steady  state  problems  because  each  of  the  iterative  solutions, 
^mrl/3^  ^nh2/3  $n+*>  satisfy  different  equations,  none  of  them  approximates 

the  steady  ctate  equations.  Even  if  the  exact  steady  state  solution  of  a 
given  problem  la  substituted  and  used  as  the  initial  data,  the  fractional 
time  step  method  will  generate  solutions  for  different  fractional  steps 
which  will  not  quite  Rettle  down  to  some  sort  of  a steady  state  Unit. 

This  situation  can  be  remedied  by  retaining  those  terms,  dropped  in  the 
fractional  step  method,  to  ba  evaluated  with  the  previous  or  otherwise 


known  iterate,  for  example. 


+ <*»  - (i  + 


At 


V 


~l  {I  + Z L 


a } (J) 


which  would  indeed  be  the  same  as  the  alternating  direction  iterative 

solution  of  the  implicit  formulation  of  the  steady  state  problem  vlth 

B ■ (I  + Lq)  end  A » — L <■  — Z Lq  given  as  equation  (4.4a) 

a [7] 

This  method  is  then  sillier  to  the  Douglas  and  Gunn's  extension  of 
Peaceman-Rachf ord' a Alternating  Direction  methods  for  the  solution  of  steady 
state  problems.  With  the  additional  terms,  it  is  not  possible  to  con- 
jecture what  the  stability  behavior  of  the  difference  formalatior  and 
the  convergence  rate  will  be.  Some  experience  of  the  research  group 
at  Langley  Research  Center,  NASA  (the  author  is  grateful  for  this  private 
communication)  Indicates  that  the  overall  computational  effort  of  using 
such  schemes  for  the  numerical  Integration  of  the  Navler-Stokes  Equations 


for  some  mixed  supersonic-subsonic  flow  fields  is  much  larger  Chen  that 
experienced  by  the  author  on  similar  problems  with  explicit  formulation. 

While  no  generalization  is  implied,  the  fundamental  reasons  expounded 
in  this  end  the  previous  sections  and  some  practical  experience  may 
serve  an  appropriate  caution  against  being  over  optimistic  to  tbe  ad- 
vantage cf  auch  Implicit  methods. 

The  fractional  step  method  was  developed  primarily  for  time  dependent 

flow  problems.  The  split  operator  L^,  for  example,  can  be  split  further  cs 

L ■ L + L where  L is  the  convective  pert  and  L the  viscous  part 
x xc  xv  sc  xv 

respectively  of  L^.  In  this  manner,  each  momentum  equation  is  split  into 
6 parts.  Each  of  the  6 parts  gives  rise  to  either  a wave  operator  or  a 
diffusion  operator.  The  question  of  rendering  a stable  computation  for  each 
step  ia  much  simplified  although  the  split  scheme  is  not  unconditionally 
stable  rnder  realistic  boundary  conditions.  There  will  be  some  inconsis- 
tencies in  the  formulation  that  can  be  remedied  by  Including  some  higher 
order  terme.  The  process  rapidly  becomes  more  complex  especially  when  high- 
er order  temporal  accuracy  is  desired!®’^  All  the  complexities  must  be 
weighed  in  the  lig^it  of  other  difficulties  in  treating  time  dependent  compu- 
tations for  3-D  flows.  Such  developments  in  dynamic  meteorology  and  ocean- 
ography can  be  important  in  aerospace  applications  in  the  near  future  if 
not  of  much  immediate  concern. 


V.  ACCURACY  AND  CONSERVATIVE  FORMULATION 


The  physical  conservation  laws  of  mass,  momentum  and  energy  are 
established  over  arbitrary  macroscopic  volumes  of  a homogeneous  fluid. 

By  reducing  the  volume  to  a macros ccpi cal ly  small  "point",  but  a micro- 
scopically large  domain  to  justify  the  continuum  model,  the  Navier-Stokes 
partial  differential  equations  were  derived  as  a convenient  mathematical 
relation  governing  smooth  point  functions  in  the  flow  field.  Now,  in  the 
interest  of  determining  some  flew  fields,  the  Navier-Stokes  equations  are 
discretized  into  a system  of  difference  equations  for  finite  elements  cf 
spatial  domains  to  facilitate  the  numerical  integration  of  the  partial 
differential  equations  system.  Such  difference  equations  may  as  well  be 
obtained  directly  from  the  consideration  of  the  fundamental  physical  laws 
for  such  finite,  discrete,  spatial  domains  with  the  help  of  interpolation 
formulas.  It  is,  however,  more  common  that  discretization  is  effected  by 
replacing  a differential  coefficient  with  difference  quotients  according 
to  some  truncated  Taylor  series  to  some  order  of  accuracy.  The  errors 
associated  with  the  interpolation  formula  or  the  truncated  Taylor  series 
ore  called  truncation  errors,  some  of  which  are  given  as  e(  in  Tables  1 
apd  II.  The  mathematical  requirement  of  consistency  means  only  that  the 
truncation  error  will  vanish  as  At,  Ax  •*•  0. 

•t 

The  conservation  laws  for  each  finite  spatial  element  arc  properly 
approximated  to  some  formal  order  of  accuracy,  according  to  the  trunca- 
tion error,  by  the  difference  equations  deduced  in  either  manner  mentioned 
above.  When  the  difference  forms  of  such  conservation  laws  are  summed 
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and  examined  over  a large  but  arbitrary  collection  of  such  finite  spatial 
elements,  however,  the  conservation  laws  may  be  seriously  violated.  This 
is  because  the  higher  order  small  errors  will  accumulate  to  much  large 
errors  when  summed  over  a very  large  number  of  the  small  discrete  elements 
which  add  up  to  the  finite  domain  of  computation.  Now  for  an  appropriate 
description  of  a physical  problem  to  the  accuracy,  say  0(Ax2),  it  is  pri- 
marily and  essentially  over  arbitrary  finite  volumes,  not  only  over  the 
differential  elements,  that  such  conservation  laws  should  b®  accurate  to 
0(Ax2).  If  the  truncation  errors  of  the  conservation  laws  in  finite  space 
is  to  be  of  0(Ax2),  it  must  not  accumulate  when  neighboring  m®sh  cells  are 
summed  up.  If  the  truncation  errors  are  allowed  to  accumulate,  the  dif- 
ference formulation  should  be  higher  order  accurate  so  that  the  accumula- 
tion of  such  higher  order  small  truncation  errors  over  arbitrary  mesh 
combinations  throughout  the  field  of  computation  will  not  exceed  0(Ax2), 
for  example.  The  difference  fora  of  Navier-Stokes  equations,  accurate 
uniformly  to  better  than  0(Ax2),  is  extremely  cumbersome  to  construct  and 
execute.  With  limited  spatial  resolution  curi'ently  available,  it  is  iter- 
ative to  prevent  or  limit  the  accumulation  of  truncation  errors. 

It  is  highly  commendable  to  verify  a posteriori,  to  what  extent  tho 
computed  results  conserve  the  mass,  momentum  and  energy  over  the  entire 
field  of  confutation.  But  this  is  not  an  alternative  to  the  requirement 
ef  no  accumulation  of  the  truncation  errors.  The  truncation  errors  are 
generally  in  the  fora  of  dipoles  or  quadruples  rather  than  siif  le  sources 
cr  sinks.  Tney  distort  the  local  flew  field  much  more  than  they  cause 
apparent  deviations  in  the  overall  mass,  momentum  and  energy  balances. 

The  consequence  of  such  dipoles  and  the  like  is  indeed  familiar  to  aero- 
dynaai cists.  A circular  cylinder  in  a uniform  incompressible  flow  is 
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represented  by  a doublet.  A thin  airfoil  or  a thin  wing  in  a subsonic 
or  a supersonic  flow  is  represented  by  some  distribution  of  sources  and 
sinks  or  dipole  pairs  within  the  framework  of  some  linearized  theory, 
known  as  the  method  of  singularities.  If  a series  of  tiny  little  vcr.es 
or  thin  sheets  are  not  to  be  tolerated  in  the  test  section  of  a windtunnel, 
such  distributed  dipoles  arising  from  the  truncation  errors  of  every  compu- 
tational cell  must  be  suppressed  if  not  completely  eliminated.  Such  sup- 
pression can  be  achieved  with  some  attention  paid  to  the  formulation  of 
the  difference  problem. 

5.1  Conservative  Difference  Formulation 

The  conservation  relations  are  written  in  divergence  fora  as  Equations 
(1.1)  to  (1.3)  for  the  density  P,  the  momentum  Pu^,  and  energy  density 
per  unit  volume.  They  are  the  scalar  components  of  the  vector  function 
V in  Equation  (1.6),  etc.  These  five  quantities  will  be  considered  as 
the  "Primary  Dependent  Variables"  in  terms  of  which  the  physical  lows  are 
stated  and  the  practical  results  desired.  They  provide  the  integrals  of 
notion  when  proper  initial  and  boundary  data  are  provided  over  a specific 
but  arbitrary  volume.  When  neighboring  volumes  are  summed,  the  contribu- 
tions on  their  common  boundary  cancel  identically  so  that  the  integrated 
conservation  laws  retain  their  identical  form.  This  is  the  crucial  pro- 
perty that  enables  the  integral  theorems  of  Stokes  and  Green  to  cast  the 
conservation  principles  into  field  descriptions  in  terms  of  different 
variables  (Fig.  2) . An  adequate  approximation  of  tho  conservation  laws 
in  difference  form  should  preferably  retain  this  property  at  least  to 
the  order  of  accuracy  required.  Such  a summable  property  is  implicit  in 
the  mathematical  abstraction  of  continuity  and  differentiability  of  the 
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functions  in  question.  Thus,  the  differential  formulations  in  terms  of 
different  dependent  and  independent  variables  are  all  equivalent  although 
the  forms  of  the  partial  differential  equations  may  be  much  different. 

This  is  not  the  case  for  the  difference  approximations  of  the  conserva- 
tion laws  that  may  be  formally  "derived”  from  the  varieties  of  forms  of 
equivalent  partial  differential  equations.  This  is  because  the  differ- 
ence functions  over  the  mesh  points  or  cells  are  discrete  or  at  least  not 
differentiable  beyond  a certain  order.  The  summable  difference  formula- 
tion, in  the  sense  that  when  cells  in  the  field  of  computation  are  summed, 
the  fluxes  in  the  physical  space  (x^)  of  the  primary  dependent  variables 

cancel  identically  along  their  common  boundary  will  be  called  "Conservative 

[3,10] 

Difference  Formulation".  The  computational  space  need  not  be  the  physical 
space.  The  dependent  variables  confuted  need  not  be  the  primary  ones. 

While  the  computation  can  be  done  in  this  manner,  it  is  still  the  fluxes 
in  the  physical  space  and  of  the  primary  variables  that  are  required  to 
be  suDDchle  for  the  conservative  difference  formulation. 

For  illustrative  purposes,  consider  tS.e  discretization  of  the  conti- 
nuity relation  from  the  integrated  conser*. aticn  law  expressed  in  the  pri- 
mary variables  p,  pu,  and  pv  in  the  two  dimensional  physical  space  (x,y) 

divided  into  miforrn  rectangular  cells  Ax  Ay.  p.  .is  the  average  density 

] ** 

of  the  fluid  in  the  cell  jAx,  kAy.  The  net  increase  of  mass  in  the  cell 
during  At  is  (p”**  - p*|  . )Ax.Ay,  The  mass  fluxes  of  pU  and  pV  should  be 
evaluated  on  the  boundary  while  pU  and  pV  are  known  only  as  the  average 
momentum  of  the  fluid  in  the  cells.  Ihus  the  boundary  fluxes  are  evaluated 
through  linear  interpolation  (but  of  second  order  accuracy)  as  the  arithe- 
matic  average  of  the  mean  momentum  in  neighboring  cells.  If  increasing 
j and  k are  positive  directions,  the  conservation  of  mass  is  stated  as: 


(5.1) 


Fcr  the  neighboring  cell  (j-l)Ax*kAy,  the  difference  form  of  mass  conti- 
nuity relation  can  be  obtained  from  (5.1)  by  replacing  j by  j-1.  The 
two  cells  have  a common  boundary  at  (j-j)Ax-kAy.  The  cut  flux  from  the 
cell  j-l,k  crossing  this  common  boundary  is  j [(pU)”  k + (pU)"_j  , which 

is  identically  the  same  as  the  in  flux  to  the  cell  (j ,k) . When  the  two 
mass  continuity  equations  (5.1)  for  the  cells  (j,k)  and  (j-1, k)  are  added, 
the  flux  terms  across  the  common  bo<mdary  cancel  out.  The  resulting 
difference  equation  is  identical  with  the  one  that  will  be  obtained  when 
the  conservation  law  is  applied  directly  to  the  combined  cells  end  is 
accurate  to  0(Ax2).  The  addition  of  other  neighboring  cells  behave  in  the 
same  manner.  Similar  results  will  be  obtained  for  the  momentum  and  the 
energy  relations.  Thus  a conservative  difference  formulation  accurate 
to  0(Ax2)  is  obtained.  It  is  easily  verified  that  the  same  difference 
formulation  will  be  obtained  with  the  forward  time,  centered  spatial 
difference  algorithm  applied  to  the  differential  equations  system  (1.1)  to 
(1.3)  written  in  divergence  form.  Indeed  the  first  order  accurate  algorithm 
of  backward  or  fcrvsrd  spatial  difference  will  also  yield  a conservation 
difference  formulation,  but  of  first  order  accuracy,  i.e.,  O(Ax)  provided 
that  tho  differential  equation  is  discretized  in  divergence  form  and  that 


the  physical  space  is  divided  into  uniform  spacing. 

If  the  continuity  equation  shbuld  be  written  in  expanded  form  for 
discretization  such  as:  u +p  17  ^or  t*ie  net  mass  ^ux  2n  x"direction, 
the  centered  space  difference  algorithm  can  represent  the  net  x-flux  as 
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Mr*  [vpm  - »m>  * pj(uj*i  - uj-i)]  <5-2a> 

or  [<UJ.l  ‘ uj.»)(pJ*l  - PJ-1>  * (pj.l  ♦ Pj-i><Vl  - uj-l)] 

(5.2b) 

Tne  in  flux  to  the  cell  (j,k)  from  the  cell  (j-l,k)  crossing  the  boundary 
at  (j-  j )Ax  is  respectively: 


^ fcjPM  + pjUj-l]  (5*3a) 

[(Vl  ♦ Uj_j)  P-.j  ♦ (pj  + 1 * P-.p  U._|]  (5.3b) 


The  out  flux  from  the  cell  (j-l,k)  into  the  cell  (j,k)  crossing  the  same 
common  boundary  as  may  be  obtained  from  Equation  (5.2a) (5.2b)  by  putting 
j -*■  j-1  is  respectively 

nr21  [uj-ipj * pj-iui]  (5-4a) 


or  JjCUj  ♦ Uj_2^Pj  * (Pj  ♦ pj-2^Uj]  (5.4b) 

The  out  flux  (5.4a)  is  identical  as  the  in  flux  (5.3a)  and  will  cancel 
each  other  when  the  two  cells  are  summed.  Thus  the  difference  algorithm 
(S.2a)  will  lead  to  a conservative  difference  formulation  with  the  differ- 
ential equation  not  written  in  divergence  form.  The  out  flux  (5.4b)  is 
different  from  the  in  flux  (5.3b).  When  the  two  cells  are  summed  up,  they 
do  not  cancel  completely  but  produce  a net  mass  source  along  the  common 
boundary,  but  in  the  interior  of  the  pair  of  cells,  with  the  magnitude 
proportional  to  Ax  . This  is  formally  negligible  in  a second  order  accurate 
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algorithm  but  renders  the  difference  formulation  from  the  algorithm  (5.2b) 
not  summable  and  not  conservative.  Even  if  such  errors  are  accumulated 
randomly  over  the  field  of  computation  with  1/Ax2  meshes,  the  accumulated 
truncation  error  will  be  of  O(Ax)  rather  than  0(Ax2).  If  the  first  order 
accurate  forward  or  backward  spatial  difference  algorithm  is  used,  the 
net  x-flux  will  be 

[UjCPj  - Pj.^  » PjtUj  - U._j)]  (5.5a) 

or  [u.Cp.^j  . P.)  . p.(a.M  - u.)]  (S.SW 


Neither  of  the  two  will  lead  to  conservative  difference  formulation  even 
at  the  accuracy  of  0(Ax).  The  above  examples  demonstrate  that  both  the 
center  difference  algorithm  and  the  divergence  form  of  the  differential 
equation  are  conducive  to  the  conservative  difference  formulation  with 
uniform  mesh  size  in  physical  space.  The  difference  formulation  based 
on  integrated  conservation  laws  even  with  linear  interpolation  on  the 
other  hand,  leads  straightforwardedly  to  conservative  difference  form 
of  second  order  accuracy. 

Consider  now  the  effect  of  mammiform  mash  3izes  in  physical  space 

with 


(Ax) 


using  the  integrated  conservation  laws  and  linear  interpolation.  The 
net  flux  into  the  cell  at  jAx  during  the  time  interval  At  is  obtained 
in  a straightforward  manner,  illustrated  here  only  for  x-fluxes: 
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At  *Ay 


nj± 

i*n 


(pu).  ♦ 


^(pu,i*V 


- AtAy 


n.  i 

±2- 


i*n 


i-i 


(pU) 


j-i 


cpw* 


(5.6) 


The  first  bracket  represents  the  out  flux  from  cell  j.  The  second  bracket 
represents  the  in  flux  to  cell  j.  If,  in  the  first  bracket,  j is  replaced 
by  j-1,  then  the  out  flux  from  the  cell  at  (Aa)^  ^ becomes  identical  as 
the  in  flux  into  the  cell  at  (Ax)^  across  their  common  boundary.  They 
cancel  each  other  when  the  two  cells  are  summed.  Thus  the  algorithm  (5.6) 
in  physical  space  will  lead  to  a conservative  difference  formulation  despite 
the  variable  spacing  in  physical  space.  The  algorithm  (5.6)  clearly  indi- 
cates how  the  centered  spatial  difference  algorithm  should  be  modified  to 
accommodate  the  variable  physical  spacing  in  order  to  achieve  the  conserv- 
ative difference  formulation  and  the  second  order  accuracy.  This  particu- 
lar combination  of  the  weighted  average  of  (pU).+j,  (pU).,  and  (pU)  . ^ is, 
however,  not  obvious  from  the  point  of  view  of  discretizing  -g-  (pu)  with 


the  second  order  accuracy  through  Taylor  series  expansions. 

It  is  common  that  variable  mesh  sizes  in  physical  space  are  accom- 
plished through  some  transf ©motion  of  the  independent  variables,  x * x(£) 
or  inversely  £ * £(x).  The  difference  formulation  is  then  derived  from 
the  transformed  differential  equation  by  discretization  over  uniform  mesh 
spacing  AC  in  the  transformed  £-space  according  to  some  difference  algorithm. 
This  transformation  of  the  spatial  coordinates  is  often  suggested  by  the 
desire  of  bringing  the  boundaries  into  coordinate  lines  such  as  £ n.  x/l+x 
so  that  x ■ • corresponds  to  £ ■ 1 or  the  use  of  spherical,  cylindrical 


or  other  convenient  body  coordinates  dictated  by  the  contour  of  the  solid 
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body  present  in  the  flow  field.  The  disretization  in  the  transformed 
C-space  in  an  intuitive  manner  is  not  likely  to  produce  a conservative 
difference  formulation.  Even  with  uniform  mesh  spacing  A£,  the  cancella- 
tion of  in  flux  and  out  flux  is  not  assured  in  the  physical  space  although 
it  is  achieved  in  the  transformed  space.  This  is  because  of  the  presence 
of  the  metric  coefficients. 

Consider  the  mass  continuity  relation  in  the  cylindrical  polar  coor- 
dinates (r,0,z): 

& * ? 4 <«“■>  * ? w <pv)  * h (pw)  ■ 0 <s-7> 

where  u,  v,  and  w are  the  radial,  azimuthal  and  axial  velocity  comncnents. 

Even  if  the  mesh  spacings  Ar,  A9,  and  Az  are  uniform  and  the  central 
space  difference  algorithm  is  adopted.  There  is  left  the  question  how 
the  metric  coefficient  r should  be  treated  in  discretizing  Equation  (5.7) 
to  obtain  a conservative  difference  form.  Now  the  integrated  conservation 
relations  in  the  physical  space  with  curvilinear  coordinates  stands  as 

Ar»Az* (r^AO) • A^ (p) 

■ At*Az  A (pur-A8)  + At*Az-ArA~(pv) 
r u 

♦ At'Ar* (r .A0) *A  (pw)  (5.8) 

J z 

where  A with  subscript  r,  0,  or  z stands  for  the  net  flux  of  the  quantity 
in  the  parenthesis  in  some  difference  form.  The  left  hand  side  of  Equation 
(5.8)  stands  for  the  net  increase  of  mass  in  the  volume  element.  If  the 
flux  terms  on  the  right  hand  side  are  expressed  either  in  the  form  (5.1) 
for  uniform  mesh  sizes  or  in  the  form  (5.6)  for  nonuniform  mesh  sizes,  for 
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example,  the  difference  form  of  the  continuity  equation  will  be  conserva- 
tive or  summable.  Thus  the  metric  coefficient  r arising  from  the  volume 
element  should  be  treated  as  r^  while  the  metric  coefficient  r arising 
from  the  surface  element  should  be  treated  differently  for  the  in  flux 
and  the  out  flux  surfaces  depending  upon  the  specific  difference  algorithm. 
It  appears  therefore  that  the  conservative  difference  formulation  can  be 
more  conveniently  obtained  by  considering  the  integrated  conservation  re- 
lations in  the  physical  space  despite  the  curvilinear  coordinate  system 
that  may  have  to  be  adopted. 

The  treatment  of  the  conservation  relations  of  the  vector  momentum 
is  considerably  more  complicated  than  that  cf  tho  scalar  mass  because  of 
the  stress  terms  and  the  inertia  terns  due  to  curvature  and  because  of 
tho  need  of  considering  tho  appropriate  vector  components . Complicated 
as  it  may  be,  the  flux  t eras  can  be  clearly  identified  and  conservative 
difference  formulations  can  be  obtained.  Often  it  is  desirable  for  achiev- 
ing a sirpler  difference  formulation  by  relaxing  the  condition  of  identi- 
cal cancellation  of  the  in  flux  and  out  flux  crossing  the  same  common  cell 
boundary.  The  more  lenient  requirement  nay  be  that  the  in  flux  and  out 
flux  crossing  the  S8ae  boundary  differ  by  a sufficiently  higher  ordar 
small  quantity  to  provide  for  some  accumulation  and  possibly  supplemented 
by  their  vanishing  over  a group  of  say  four  neighboring  cells.  This  may 
be  permissible  since  the  ultimate  objective  of  the  conservative  differ- 
ence formulation  is  to  prevent  undue  accumulation  of  truncation  errors 
over  finite  volumes  to  cause  serious  deterioration  of  the  accuracy  of  the 
confutation. 


>WWir; 
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With  conservative  difference  formulation,  the  accumulated  truncation 
error  Ej.  of  a set  of  calculations  can  be  estimated  to  the  order  of  magni- 
tude at  any  point  within  the  field  of  computation.  Moreover,  the  error 
of  the  computed  results  at  a point  can  tc  separated  into  two  parts: 

(i)  the  tnmeatien  error  Ey,  and  (ii)  the  boundary  error  at  the  point 
caused  by  the  errors  on  the  boundary  of  the  field  of  computation  although 
the  difference  problem  is  essentially  non-linear. 

S.2  Heuristic  Error  Estimate  and  Accuracy 

For  a computed  result  to  be  practically  useful  it  is  essential  to 
have  some  idea  of  its  accuracy.  The  study  of  the  accuracy  question  has, 
however,  been  little  explored.  This  may  be  due  partly  to  the  preoccupa- 
tion with  stability  questions  and  partly  to  the  difficulty  of  construct- 
ing an  upper  bound  of  the  error  of  a computation  for  the  type  of  initial 
boundary  value  problems  of  fluid  dynamics.  It  may  be  possible  that  those 
convergence  proofs  which  naturally  include  the  estimate  of  the  error 
bounds  will  be  extended  from  the  periodic  boundary  value  problems  to  more 
realistic  boundary  conditions  as  applied  to  the  Navier  Stokes  equations. 

In  practice  such  a difficult  and  complicated  a-priori  rigorous  error  esti- 
mate is  not  necessary.  Heuristic,  a-posteriori  rough  error  estimates 
will  often  suffice.  Indeed,  it  would  be  preferable  to  have  the  estimate 
simple  and  generally  applicable  although  not  rigorous  and  not  so  precise. 
Tho  nonlinear  Burgers'  equation  is  therefore  conveniently  adopted  for 
analysis  as  a one-dimensional  model  of  the  Navier-Stokes  equations! 10 Us 
was  shewn  in  Section  (3.2),  it  is  a useful  model  for  stability  analysis, 
being  quasi- linear  and  with  both  wave  and  diffusion  characteristics.  It 
is  also  convenient  for  the  study  of  accuracy  becauso  many  exact  solutions 


78  - 


are  known,  with  which  the  coinputational  errors  can  be  quantitatively  eval- 
uated and  compared  with  theoretical  estimates. 

The  Burgers'  equation,  in  dimensionless  form  is 


3u 

7t 


♦ u 


3u 

cix 


1 32u 


with  a steady  state  solution 


(5.9) 


u(x)  - - a tail  (aRex/2) 


(5.10) 


having  u(x*0)  • 0 u(x=-l/2)  ■ 1 


and  |u(x»±»)|  » a = 1/tanh  (aR  /4) 

c 

This  steady  state  solution  in  the  range  of  -1/2  <.  x <.0  has  been  calculated 
as  the  long  time  limit  of  the  temporal  problem  via  several  difference 
algorithms.  Hie  quasi-linear  term  u^^-  is  always  treated  in  the  divergence 
fora  (u2/2)  with 


U1  /2(2-*-a) 


(5.11) 


and 


Here  "a"  is  a parameter.  The  single  centered  spatial  difference  corres- 
ponds to  a - 0.  The  center  spatial  difference  in  non -divergence  form 
results  when  a « «,  in  which  case. 


’ ujLL 


2Ax 


) 


(5.12) 
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If  it  is  presumed  that  an  approximate  steady  state  solution 
U(t,x,At,Ax)  will  be  reached,  which  departs  only  slightly  from  the  gen- 
uine solution(5.10) , a linearized  differential  equation  for  the  error 
can  be  derived  and  solved.  The  linearization  permits  the  separation 
of  the  errors  as  originated  from  various  sources,  i.e.,  as  the  sum  of 
the  truncation  errors  and  the  boundary  errors  E^  The  difference 
equations  derived  from  (5.11)  are  all  nonlinear  but  conservative  and 

permit  the  separation  of  the  truncation  and  the  boundary  errors  with 

/ 

the  cumulative  truncation  errors  ET  remain  of  the  order  of  tho  local 
truncation  error. 

The  linearized  snalysis  shews  that  Ep  at  any  point  in  the  field 
of  ccs$>utation  is  proportional  tofRe^)2  for  the  second  order  accurate 
conservative  difference  formulations  derived  from  (5.11).  Here  Re^x  is 
the  Reynolds  number  based  on  the  length  Ax  and  the  velocity  difference 
between  the  point  x-0  with  maximum  velocity  gradient  end  the  point  x«l 
noarly  tho  asymptotic  velocity.  For  a quantitative  estimate  of  the  E_, 


ME.  ♦ 


" "ofco 


MjEj  ♦ (l*3a)M2E2  Mj 


2 + a 


♦rEs 


(5.13) 


where  MQ  is  the  constant  defining  the  steady  state  criterion 


(R®£X) 2 and  M^^Ax^2  Rre  th®  co«fficicnts  of  the  truncated  quasi- 
linesr  convective  terms  end  M3(Re^x)2  is  that  of  the  truncated  viscous 
teras.  Mj,  M2,  and  Mj  are  expected  to  be  of  0(1)  for  reasonable  differ- 
ent algorithms  and  for  reasonably  smooth  solutions.  E , E.,  and  E_  are 
universal  functions  of  the  genuine  solution  u(x)  that  vanish  on  both 
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bound3ri.es  and  h aye  their  absolute  magnitudes  less  than  0.1.  Thus  the 
truncation  errors  E^  are  expected  to  be  of  the  order  of  (ReAx)2/10  for 
the  second  order  accurate  schemes.  Actual  computations  with  Mq  ■ OCAx)1* 
and  Ax  « 1/20  for  various  schemes  verified  tho  quantitative  values  of 
Equation  (5.13)  end  the  dependence  of  on  CRe^x)** 

For  Re^x  ■ 0(1)  and  all  the  finite  values  of  a * 0(1)  tested,  the 
following  estimate  of  the  maximum  absolute  truncation  errors  is  valid 

Et  < 3 x 10" 2 (Re^)*  (5.14) 

This  simple  formula  is  therefore  recoaanended  es  a preliminary  estimate 
of  the  bounds  of  the  truncation  errors  of  second  order  accurate  conserv- 
ative difference  formulation.  With  non-conservative  difference  formula- 
tion, the  truncation  errors  can  accumulate  and  become  considerably  larger 
than  the  estimate  given  by  Equation  (5.14). 

The  boundary  errors  in  tho  field  due  to  a fractional  error  of  the 
boundary  value  is  given  by  the  linearized  analysis  as: 

ED-ebEh  (5.15) 

where  is  t universal  function  that  is  unity  on  the  boundary  where  the 
erroneous  boundary  condition  is  applied,  end  decays  very  slowly  toward 
the  other  boundary  where  it  vanishes.  The  decay  is  so  slow  that  the  error 
retains  more  than  half  its  value  until  within  the  last  tenths  of  the 
field  of  computation  near  the  other  boundary  (note  that  is  plotted  against 
u(x)  in  Fig.  2 of  Ref.  10)  depending  upen  the  magnitude  of  the  Reynolds  num- 
ber. 

For  Neumann  boundary  conditions,  the  boundary  error  is  still  given  by 
Equation  (5.15)  but  is  evaluated  as 
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Ej,  - -2c^/aRe  (5.16) 

where  is  the  fractional  error  of  the  spatial  derivative  on  the  boundary. 
Within  the  framework  of  linearized  estimate  of  errors,  the  superposition 
of  (5.15)  and  (5.16)  with  proper  coefficients  will  enable  an  estimate  of 
the  errors  caused  by  a Cauchy- type  condition.  The  boundary  error  at  a 
given  point  in  the  field  of  computation  will  be  the  sum  of  the  decayed 
boundary  errors  from  both  boundaries. 

In  multidimensional  flow  problems,  it  is  presumed  that  the  results 
of  the  previous  model  analysis  may  apply  primarily  in  the  direction  along 
streamlines  or  nearly  so.  This  leaves  the  estimate  of  the  contributions 
of  the  boundary  errors,  from  those  portions  of  the  boundary  of  the  field 
of  computation  that  are  primarily  parallel  to  the  local  streamline  direc- 
tions, yet  to  be  accommodated.  No  helpful  suggestions  can  be  made  here 
except  to  render  a description  as  nearly  correct  as  the  physical  situation 
suggests.  In  fact,  the  treatment  of  this  portion  of  the  computational 
boundary  is  one  of  the  two  outstanding  difficulties  that  the  author  and 
his  co-worksrs  have  experienced  in  various  problems.  (The  other  outstand- 
ing difficulty  is  the  treatment  of  internal  shockwaves  to  be  explored  in 
tho  next  section.) 

The  decay  characteristics  described  by  the  universal  function  F^ 
may  be  used  where  the  one -dimensional  model  is  appropriate.  The  various 
universal  functions  Eq,  Ej,  etc.  and  E^  in  the  model  results  may  bo 
recognized  as  the  "influence  functions"  describing  the  error  propagation 
in  the  field  of  computation.  They  can  be  empirically  established,  a 
posteriori,  by  introducing  a known  error  at  a specific  point  (on  the 
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for  truncation  errors)  and  then  computing  the  solution  under  the  modulated 
condition.  The  difference  of  the  two  sets  of  solutions  then  gives  the 
influsnce  function  E in  question.  Usually,  during  the  development  stage 
of  a difference  formulation  for  a given  physical  flow  problem,  such  infor- 
mations can  be  derived  from  preliminary  results  and  can  be  used  for  the 
purpose  of  a-posteriori  error  estimate.  Of  course,  a-posteriori  determi- 
nating of  such  influence  functions  are  desirable  to  provide  additional 
checks  on  the  behavior  of  the  computational  program. 

Without  referring  to  any  specific  computational  problem,  the  follow- 
ing general  observations  can  be  inferred  from  the  model  study.  They  are 
applicable  only  for  the  conservative  difference  formulation  in  which  the 
truncation  errors  do  not  accumulate  so  that  the  truncation  and  the  boundary 
errors  can  be  treated  separately  and  estimated  by  Equations  (5.14)  to  (S.16) 

(1)  The  steady  state  criterion  |u*?+1  - U^l  < O(Ax)1*  is  sufficiently 
accurate  in  a second  order  accurate  scheme. 

(2)  The  truncation  error  E^,  is  expected  to  be  MRe^x)n  ^or  conserva" 
tive  difference  formulation  of  nth  order  formal  accuracy  and  the  influence 
fractions  Ej  etc*  aye  not  likely  to  possess  maximum  magnitudes  much 
less  than  10”*.  With  Re^x  > 1 in  the  practical  cases,  maximum  truncation 
error  is  not  likely  reduced  appreciably  from  that  of  second  order  accurate 
scheme  as  may  be  estimated  from  (5.14). 

(3)  Boundary  errors  cannot  be  efficiently  reduced  by  reducing  the  mesh 
sixes.  They  decay  vary  slowly  and  are  generally  considerably  larger  than 
the  truncation  errors  in  practical  cases  with  Re^x  » 0(1).  The  primary 
effort  required  in  achieving  a reasonably  accurate  solution  of  the  conpli- 
cated  practical  problems  lies  in  the  sophistication  in  the  treatment  of 

the  various  boundary  conditions.  The  field  of  computation  and  the  choice 
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of  coordinates  should  be  properly  defined  to  facilitate  a Bore  accurate 
iop  lamentation  of  the  boundary  conditions. 

Hie  general  observations  made  above  carry  an  important  message  to 
those  interested  in  obtaining  solutions  for  the  complicated  fluid  dynam- 
ics problems  with  reasonable  accuracy  to  suit  practical  purposes.  Much 
attention  should  be  paid  to  the  formulation  of  the  problem.  Attempts  to 
improve  the  accuracy  of  a numerical  solution  based  on  a poor  formulation 
by  extending  the  computation  to  satisfy  more  restrictive  steady  state 
criterion,  with  more  refined  mesh  and  even  with  the  help  of  much  larger 
and  faster  coaputers  can  prove  to  be  not  only  expensi’/e  but  frustrating. 

A similar  attempt  was  made  for  the  time  dependent  flow.  It  was  found 
that  for  flows  with  slow  and  monotonic  temporal  variations,  the  behavior 
of  error  propagation  in  the  sscond  order  accurate  conservative  difference 
formulation  is  essentially  similar  to  that  described  above  for  steady 
state  problems.  For  oscillatory  flows,  conservation  in  the  spatial  space 
apparently  fails  to  holp.  Test  calculations  for  some  simple  damped  oscil- 
lation as  an  exact  solution  of  the  Burgers'  equation  indicate  the  sorious 
effect  of  the  phase  errors  of  different  oscillatory  components  caused  by 
the  dispersive  truncated  terms.  It  has  been  suggested  and  illustrated 
that  fourth  order  accurate  difference  algorithms  will  substantially  inprovc 
the  accuracy  of  the  computed  results  beyond  quite  a few  cycles  of  oscil- 
lations. It  is,  however,  a tremendous  task  to  compute  with  fourth  order 
formal  accuracy  and  uniformly  for  as  complicated  an  equation  system  as  the 
Navier-S  tokes J 8 »9  * 11 1 
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5.3  Shock  Wave  and  Artificial  Viscosity 

In  all  the  previous  discussions,  the  question  of  "non-smooth"  or 
even  "discontinuous"  solutions  is  deferred.  In  problems  of  practical 
interest,  shock  waves  and  contact  discontinuities  are  often  the  promi- 
nent features  of  the  flow  problems.  The  presence  of  such  discontinuities, 
or  generally  regions  of  very  large  gradient,  causes  difficulties  in  their 
coaputation. 

Discontinuous  initial  and  boundary  data  ere  often  imposed  on  purely 
elliptic  or  parabolic  problems.  They  may  cause  oscillations  in  the 
vicinity  of  the  boundary,  but  never  very  serious.  This  is  because  of  the 
inherent  nature  of  the  system  to  smooth  out  any  discontinuities  in  time 
and  in  space.  The  accuracy  of  the  computed  results  may  suffer  somewhat 
according  to  the  modulus  of  continuity  of  the  functions  involved,  but  can 
often  be  remedied  by  using  a higher  order  accurate  difference  algorithm. 

This  inherent  tendency  to  smooth  out  any  discontinuity  can  also  be  trouble- 
some, for  example,  in  treating  flow  problems  involving  an  interfacial 
discontinuity  formed  by  two  different  fluid  mediums,  especially  when  the 
interface  is  not  stationary;  since  initially  sharp  discontinuity  dif- 
fuses in  the  course  of  the  coaputation  if  not  artificially  maintained. 

For  hyperbolic  problems,  a discontinuity  in  the  initial  bowdary 
data  propagetes  into  the  field  of  computation  and  causes  excessive  compu- 
tational disturbances  downstream,  particularly  in  its  zone  of  influence. 

It  also  produces  upstream  influences.  For  the  quasi- linear  gas  dynamic 
problems,  a shock  discontinuity  can  physically  arise  from  a perfectly 
smooth  boundary  due  to  the  coalesence  of  the  smooth  compression  waves. 

Thus,  when  this  flow  field  is  computed  with  an  algorithm  that  worked  well 
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for  smooth  fields,  quite  severe  oscillations  can  develop  at  approximately 
where  the  shock  discontinuity  would  appear.  Such  oscillations  are  fairly 
large,  although  not  necessarily  leading  to  the  catastrophic  divergence  of 
linear  instability.  It  may  be  that  the  amplitude  of  such  shock  induced 
oscillations  are  limited  by  the  non-linear  effects  and  the  phenomenon  may 
well  be  called  nonlinear-instability,  but  certainly  not  in  the  sense  of 
violating  the  requirement  of  boundedness  discussed  in  Chapters  II  and  III. 
Even  if  boicided,  such  oscillations  are  highly  damaging  to  the  accuracy  of 
the  results  not  only  in  the  vicinity  of  the  shock  but  over  most  of  the  flow 
field  computed.  Practical  interest,  on  the  other  hand,  is  often  centered 
in  the  vicinity  of  such  ahock  discontinuity. 

It  is  natural  to  treat  the  shock  front  or  the  interfacial  contact 
discontinuity  as  an  internal  boundary  and  to  compute  the  smooth  solutions 
on  both  sides  of  the  discontinuity  separately.  The  jump  conditions  across 
the  discontinuity  will  connect  the  two  solutions  together.  This  shock 
matching  or  shock  fitting  procedure  is  easily  carried  out  in  one  space 
dimension  for  a known  discontinuity,  i.e.,  the  magnitudes  of  the  jumps 
and  the  speed  of  propagation  of  the  discontinuous  front  into  a homogene- 
ous medium  at  rest  or  in  uniform  motion.  If  the  shock  should  be  propagat- 
ing into  a non-iaiifora  medium  or  a homogeneous  medium  in  non-uniform 
motion,  the  shock  strength  and  speed  will  vary  and  the  Hugoniot  relations 
across  the  shock  will  have  to  be  s implemented  by  some  additional  matching 
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vicinity.  The  oscillations  may  be  alleviated  if  the  shock  location  is 
fixed  at  a mesh  point  and  the  mesh  divisions  are  rezoned  at  every  time 
or  iterative  steps.  The  computational  procedure  in  terms  of  such  shock 
coordinate  rapidly  becomes  complicated. 

In  two  space  dimensions  and  with  a curved  shock  of  unknown  shape  and 
location,  the  computational  details  of  such  a shock  matching  procedure 
rapidly  becomes  more  tedious  and  inaccurate.  With  fixed  mesh  points,  the 
shock  front  is  generally  off  the  mesh  points  and  it  becomes  difficult  and 
highly  inaccurate  to  determine  the  direction  normal  to  the  front  in  the 
matching  process.  The  use  of  curvilinear  shock  coordinate  is  convenient 
and  may  possess  other  features  for  treating  inviscid  steady  state  flow 
problems  with  uniform  suporsonic  flow  on  the  upstream  side  of  the  shock 
feont<.12lt  is  not  suitable,  however,  for  a shock  wave  ishedded  in  a non- 
uniform  inviscid  flow  field.  It  cannot  be  implemented  for  viscous  and 
inviscid  flow  fields  involving  more  complicated  shock  configurations, 
such  93  shock  intersections  and  Mach  reflections  or  transonic  shocks  that 
terminate  in  the  flew  field.  The  tedious  shock  matching  can  be  imple- 
mented in  principle  even  for  such  complicated  configurations  although  the 
procedure  is  too  complicated  to  be  manageable  and  the  results  so  obtained 
are  uniformly  poor. 

[13] 

To  avoid  shock  matching,  v. Neumann  and  Richtoeyer  introduced  the 
artificial  viscosity  uothod  for  cosqputing  the  shock  propagation  in  an 
inviscid  flow  field.  A quadratic  viscous  pressure  term  pa2Ax2 , 
where  a is  a numerical  constant  chosen  at  convenience,  is  added  to  the 
differential  equation  before  discretization.  The  quadratic  dependence 
on  the  velocity  gradient  assures  a rapid  decay  of  the  artificial  viscous 
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term  away  from  the  shock  front  with  steep  velocity  gradient.  With  a < 1 
the  typical  results  of  the  calculation  for  one-dimensional  shock  propa- 
gation into  a uniform  field  gives  a sharp  shock  front,  spreading  over 
n.2  meshes  and  the  calculated  shock  speed  is  within  0.1%  of  the  correct 
value.  But  sizable  oscillations  develop  in  the  downstream  and  over  an 
extended  range  without  appreciable  damping  (spatially  end  teaporally) . 

By  increasing  a to  £ 2,  the  magnitudes  of  the  oscillations  are  reduced 
but  the  shock  front  spreads  wider,  over  4 or  more  meshes.  A reasonably 
smooth  downstreem  solution  is  obtained  only  when  a is  so  larjp  as  to  be 
0(Ax_1)  and  the  shock  front  spreads  over  many  mashes.  By  then  the  arti- 
ficial viscous  terra  is  no  longer  smell  in  the  inviscid  region  and  the 
apparently  smooth  results  of  computation  fail  to  be  a satisfactory  approx- 
imate solution  for  the  shock  front. 

The  artificial  viscosity  method  is  physically  sound,  simply  imple- 
mented and  easily  extended  to  multispace  dimensions,  formally  by  includ- 
ing dorivatives  in  other  spatial  dimensions.  The  large  spread  of  the 
shock  front  and  the  induced  oscillations  generally  become  more  objection- 
able, however.  Many  artifices  can  and  have  been  devised  to  improve  the 
appearance  of  the  computed  results.  The  artificial  viscous  term  may  be 
dropped  when  the  gradient  of  velocity  becomes  less  than  a pre-assigned 
value,  or  the  downstream  oscillations  may  be  suppressed  or  eliminated  by 
some  smoothing  process  or  may  be  limited  to  a permissible  magnitude  about 
the  mean  through  some  filtering  process.  Excellent  results  can  generally 
be  obtained  for  simple  test  problems  with  known  shocks.  The  merit  of 
such  procedures  in  computing  shock  propagation  into  non-uiifora  flow 
fields  is  yet  to  be  demonstrated,  particularly  with  respect  to  the  accui  zy 
of  the  smooth  results. 
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The  Lax-Wendroff  treatment  of  shock  wave  utilizes  the  fact  that  the 
hugoniot  relations  are  simply  the  conservation  laws  integrated  over  the 
discontinuity.  Thus,  with  the  inviscid  equations  written  in  divergence 
fora  for  the  physically  conserved  quantities,  shock  matching  can  be 
avoided  since  the  difference  equations  for  such  conserved  quantities  are 
indeed  the  approximate  form  of  the  hugoniot  relations.  (Note  that  the 
divergence  form  of  the  transformed  dependent  variables  does  not  help.) 

One  dimensional  computations  show  that  it  leads  to  a quite  sharp  shock 
front  (y2bx)  and  accurate  shock  speed.  But  sizable  oscillations  are 
generated  at  the  shock  frent  although  rapidly  daa^ped  to  disappear  within 
8 to  10  meshes  from  the  front.  Tho  dissipation  is  derived  from  the  dissi 
pative  term  Ax3»r(l-r2)i£j^-  with  r * uAt/Ax  which  may  be  visualized  as  the 
artificial  viscosity  that  spreads  out  the  shock  front.  The  quadratic 
viscou3  terms  adopted  by  v. Neumann  and  Richtmeyer  does  not  appear  to  pro- 
vide as  much  damping  of  the  shock  induced  oscillations  as  this  linear  vis 
coue  term.  The  peak  ajqjlitude  of  this  shock  induced  oscillation  near  the 
front  is  often  larger  although  more  rapidly  damped  than  those  from  the 
quadratic  artificial  viscous  term.  Additional  artificial  viscous  terms 
aro  often  introduced  to  reduce  the  amplitude  of  such  oscillations. 

The  introduction  of  artificial  viscous  terms  into  the  differential 
equation  before  discretization  is  fundamentally  not  much  different  from 
dropping  those  higher  order  terms  in  the  truncated  Taylor  series  in  the 
discretization  process.  Such  viscous  terms  contribute  to  the  stability 
of  the  difference  formulation.  Thus  artificial  viscosity  is  very  widely 
employed  for  problems  without  shocks.  Such  artificially  introduced 
viscous  terms  are  often  substantially  larger  than  the  Navier-Stokes 
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viscous  stress  terns  evaluated  with  the  physical  viscosity  coefficient 
of  the  fluid.  This  is  justifiable  for  the  solution  of  inviscid  flow 
problems,  i.e.,  the  flow  problems  visualized  as  the  asymptotic  lirdt  of 
negligible  viscous  stress  terms,  so  long  as  the  contributions  of  the 
artificial  viscous  terms  are  "negligibly  small"  compared  with  the  contri- 
butions from  the  inviscid  terms  and  the  somewhat  spread  out  shock  front 
is  visualized  as  a "sharp"  discontinuity.  Such  a situation  is  clearly 
not  tolerable  for  viscous  flow  problems  because  the  effect  of  the  fluid 
viscosity  will  be  overshadowed  by  the  effect  of  the  pseudo  viscosity. 

There  are  many  numerical  solutions  of  the  Navier-Stokcs  equations, 
some  with  first  order  accurate  algorithms,  some  with  second  order  accur- 
ate algorithms  but  with  large  artificial  viscous  terms,  at  largo  Reynolds 
nunfcers  based  upon  fluid  viscosity  of  the  order  of  10^.  These  computed 
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results  are  very  insensitive  to  the  large  fluid  Reynolds  numbers.  This 
is  understandable  since  the  pseudo- viscosity  in  Such  calculations  are 
substantially  larger  than  the  real  fluid  viscosity,  and  changes  in  the 
fluid  Reynolds  number  will  not  significantly  alter  the  effective  Reynolds 
number  of  the  conputational  results  based  on  the  total  viscosity  included 
in  the  difference  formulation.  If  one  wishes  to  evaluate  quantitatively 
the  viscous  effects,  both  the  artificial  viscous  term  introduced  into  the 
differential  equations  system  and  the  pseudo  viscous  terms  implicit  in  the 
difference  form  should  remain  substantially  less  than  the  fluid  viscous 
term.  Thus  for  viscous  flow  problems,  artificial  viscosity  terms  of  the 
type  used  by  v. Neumann  and  Richtmcyer  should  satisfy 
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or  dimensionally 


» AuAx 
v 


Re4x  I 


(5.17) 


With  Re^  generally  larger  than  unity,  the  constant  a must  be  chosen 
appreciably  less  than  unity.  This  restricts  severely  the  usefulness  of 
the  idea  of  using  the  artificial  viscous  term  either  for  securing  compu- 
tational stability  or  for  suppressing  the  shock  induced  oscillations  in 
viscous  flow  problems. 

For  second  order  accurate  conservative  difference  formulation,  the 
errors  introduced  by  the  pseudo- viscous  terms  are  included  in  the  trunca- 
tion error  ET,  the  absolute  upper  bound  of  which  may  be  estimated  as 
< 3 x 10  ^ReAx^2  accor<^n2  t0  fhe  results  based  on  the  Burgers' 
model  equation  given  in  the  previous  section.  Thus  Re^  may  be  as  large 
as  1 or  even  2 without  having  the  cumulative  truncation  errors  exceeding 
a few  percent.  Note  that  this  Re^x  is  defined  in  terms  of  the  local  change 
cf  volocity  per  mesh  when  the  Burgers'  model  is  applied  to  the  "local 

flow  field".  With  the  Reynolds  number  based  on  the  viscous  flow  dimension 

3 4 2 

and  maximum  velocity  of  0(10  - 10  ) , and  with  possibly  of  0(10  ) mesh 

points  over  the  linear  dimension,  the  local  values  of  Re^  will  generally 

be  considerably  smaller  than  10  (since  Au  per  mesh  will  be  significantly 

less  than  10  *)  except  possibly  in  the  region  of  shock  induced  oscillations. 

If  the  shock  front  is  'isualized  as  an  interior  boundary  and  the  shcck-in- 

duccd  oscillation  as  a form  of  propagating  boundary  error,  the  errors  in 

the  results  computed  with  the  second  order  accurate  difference  formulation 

will  generally  be  dominated  by  boundary  errors. 
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Shock  induced  oscillations  mar  the  appearance  of  the  computed  solu- 
tion nuch  nore  seriously  although  need  not  cause  larger  error  than  the 
less  conspicuous  sources  froa  the  exterior  boundary.  The  difficulty  is 
compounded  where  a shock  wave,  either  incident  or  emerging,  intersects 
the  exterior  boundary.  In  the  next  section,  the  relation  between  the 
boundary  treataent  and  the  shock-induced  oscillation  will  be  explored. 

5.4  Shock-Induced  Oscillations 

Shock  induced  oscillations  are  often  considered  as  unavoidable  when 
a shock  wave  is  encowitered  in  the  confutation  with  a higher  erdsr  accur- 
ate difference  algoritha.  While  the  first  order  accurate  algorithm  will 
not  give  rise  to  such  oscillations,  the  smear  of  the  shock  front  becomes 
excessive  and  the  emulative  truncation  errors  becomes  large.  Thus  when 
shock  wave  is  encountered  in  a confutation,  it  is  often  held  as  necessary 
to  choose  between  the  two  evils.  The  following  is  an  attempt  to  clarify 
the  origin  of  the  spurious  oscillations  and  to  show  that  a certain  class 
of  second  order  accurate  difference  algorithms  can,  under  some  favorable 
circumstances,  avoid  such  spurious  shock-induced  oscillations. 

Consider  the  solution  of  a linear  steady  state  problem  via  the  time 
dependent  approach.  Let  the  spatial  difference  operator  be  split  into 
. two  parts,  L^CT)  and  L2(T)  where  T is  the  shift  operator  for  the  spatial 
indices,  i.e.,  TU^  - Uj+1,  T"Vj  - U^j,  and  tV  - T-Tlh  « Uj+2,  etc. 
Construct  the  class  of  two  step  difference  algorithms  for  the  time  inter- 
val nAt  to  (n-tl)At: 

Uj  - . LjCnu*?  ♦ L2(T)uJ 

Uj+1  - - LjCT)^  ♦ L2(T)UJ  (5.18) 
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where  if?  is  a provisional  or  predicated  value  of  U*?**  and  the  second  or 
final  step  is  a corrector  step.  LjfT)  ♦ L2(T)  is  second  order  accurate 
and  consistent  with  the  differential  operator  in  the  steady  state. 

Let  the  boundary  conditions  to  be  appliod  to  the  first  or  the  pro- 
visional step  be 

B(T)  ty?  - 0 (5.19) 

and  let  the  boundary  values  of  if?, derived  from  these  boundary  conditions 
as  are  used  in  the  first  stejv  be  used  in  the  second  step  for  the  computa- 
tion of  at  the  corresponding  boundary  points.  In  this  manner  it  is 
maintained  that  • 1j?  = 0 on  all  the  boundary  points  and  at  every 
tine  step.  The  boundary  values  at  each  boundary  point  may  change  from 
step  to  step  and  contain  errors  implicit  in  the  boundary  conditions  (5.19). 

By  subtracting  the  two  steps  in  the  difference  equations  (5.18)  the  follow- 
ing difference  relation  is  obtained: 


In  the  event  that  a steady  state  is  approached  in  the  sense  that 
Uj*1-  uj.  then  Equation  (5.20)  becomes  in  the  steady  state  limit 


[«♦  v»]f5-'5) 


(5.20) 


(S.21) 


Thus  - U**  is  governed  by  the  linear  system  of  difference  equations 
(5.21)  and  are  subject  to  zero  boundary  values  over  the  entire  boundary. 

If  there  is  no  eigen  solutions  to  this  system  of  equations,  it  follows  in 
the  steady  state  limit  that  Then  the  solution  in  the  steady 
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state  limit  is  the  solution  of  the  correct  steady  state  equation 

[l^t)  ♦ l2(t)]  - 0 

Now  if  the  boundary  values  of  ^f?  and  are  not  kept  the  same  in 
successive  iterations,  if?  will  have  to  be  eliminated  from  Equations  (5.18). 
That, in  the  limit  of  the  steady  state  with  ifj*1  » U*j,  this  solution  will  be 
determined  by  the  equation 


(5.22) 


£l  ♦ Lj(T)  ]f  j(T)  ♦ L2(T) J Uj  - 0 


(5.23) 


It  will  contain  the  "correct"  steady  state  solution  (5.22)  to  the  extent 
that  the  boundary  conditions  B(T)lf)  - 0 represent  the  correctly  posed 
situation.  But  it  also  contains  the  nontrivial  solutions  of  Equation 
(5.21)  when  ifj  - if?  is  not  identically  zero,  as  a result  of  the  slight 
difference  in  the  boundary  values  of  and  Such  extrane- 

ous solutions  naturally  are  possible  sources  of  the  shock-induced  oscil- 
lations end  can  indeed  be  identified  in  the  course  of  co amputation  as 
being  proportional  to  the  difference  between  the  provisional  and  the 
final  solutions.  From  the  practical  point  of  view,  it  is  simpler  and  more 
desirable  to  use  identical  boundary  values  from  (5.19)  to  suppress  all  the 
spurious  fuidaamtal  solutions  arising  from  Equation  (5.21). 

There  are  many  two  step  difference  algorithms,  but  mostly  not  of  the 
class  (S.18)  except  the  Cheng- Allen  scheme  and  the  Brailovskaya's 
schema.  For  the  linearized  Burgers'  Equation  (3.13),  the  difference 
fores  can  be  cast  into:  Cheng -Allen  Algorithm^0' 16 ^ 


•rar[(-T‘5)T*(r*,)T‘1] 


l2(t) 


-2s 

TTZs" 


(5.24) 
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Brailovskaya  Algorithm 


U7] 


Lj(T)  - r(T  - T'1) 
L,(T)  - s(T  - 2 ♦ T'1) 


where  r ■ cAt/Ax  and  s • vAt /Ax2.  When  (5.24)  is  substituted  into 
Equation  (5.23),  the  general  solution  of  is  obtained  as 


(5.25) 


where 


and 


Uj  mL\  k*  k " 1.  2.  3.  < 


^l-1 


2*  »T 

is-r 


**7  P°Ax 

lT“*Ax 


53  4 - £(!♦*).*  j(l*2s)*  ♦ (ra-sl)}l/2]^ 


(2s-r) 


(5.26) 


and  are  the  two  proper  fundacentnl  solutions  of  the  correct  steady 
state  equation  Q-jCT)  ♦ L2(T)JUj  ■ 0,  because  in  the  limit  of  Re^x  -*■  0, 
they  approach  the  two  fundamental  solutions  1 and  exp (Rex)  of  the  steady 
state  differential  equation,  c|j  ■ ]£■ -jjr . £3*  an<*  £4  ,re  two  extran- 
eous fundamental  solutions  of  the  two-step  scheme  that  constitute  the 
errors  or  "spurious  solutions"  arising  from  the  solution  of  the  equation 


E * iim]uj 


or  of  Equation  (5.21). 
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Kith  both  r and  s > 0 , and  ‘^y1-  < 1 , it  is  found  that 


*4 < 0 “ r i *• 


Thus  5^  always  represents  mesh  to  mesh  oscillation  while  ^ can  be  either 
oscillatory  or  monotomic.  Hie  steady  state  licdt  of  the  difference 
if1*1  can  be  given  as: 


V V W * C454 


j 


(5.27) 


where  Cj  and  c^  are  deternined  by  the  difference  in  the  boundary  values 
of  > and  iT1  at  j ■ 0 and  j « J cn  the  boundary.  When  the  boundary 
values  of  u”  and  U***1  are  kept  the  same  at  every  step,  then  c^  ■ c^  ■ 0 
and  no  spurious  solution  will  be  present  in  the  computed  steady  state 
result.  Otherwise,  oscillations  con  be  expected  in  the  computed  steady 
state  result. 

If  Brailovskaya's  scheme  (5.25)  is  substituted  into  Equation  (5.23), 
the  sane  proper  fundamental  solutions  Sj-*  and  ^ tre  obtained,  but  the 
pair  of  extraneous  solutions  and  ^ are  given  somewhat  differently  ' 
with  ^ » [l  i (1  ♦ 4r2)1^2  J/2r,  and  < 0 always.  The  overall  situ- 
ation is  much  the  same. 

It  may  be  pertinent  to  repeat  here  that  the  spurious  solutions  will 
be  suppressed  so  long  as  the  same  boundary  values  of  ?f*  and  u”*1  are  used 
at  every  step.  Such  boundary  values  can  be  detersined  by  the  approximate 


boundary  conditions  B(T)U^  » 0,  and  may  contain  errors.  In  this  event. 


they  nay  cause  errors  in  the  constants  Cj  and  c2  of  the  steady  state 
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solution 


U,  - c1?1J  . c2?2i 


(5.28) 


There  will  not  be  any  catastrophe  if  the  boundary  values  are  not  exces- 
sively in  error  and  if  the  mash  size  of  the  steady  state  solution  is 
not  too  coarse  so  that 


R*Ax  < 2 (5.2S 

This  last  restriction  Re^x  < 2 has  little  to  do  with  suppressing 
the  spurious  fundamental  solution,  m<1  It  is  to  keep  from 

being  oscillatory  and  failing  to  be  a valid  approximation  to  the  fundamental 
solution  exp  |RejAx|  of  the  differential  problem.  It  is  clear  from  Equation 
(5.26)  that  when  Re^x  > 2,  the  appropriate  form  of  ^ is 


t-”3  |U 


which  is  oscillatory  and  rapidly  amplifying  with  increasing  j.  It  fails 
to  serve  as  any  meaningful  approximation  to  exp  |RejAx|.  Thus  to  obtain 
a valid  steady  state  solution  without  spurious  oscillations  based  on  the 
algorithms  (5.24)  or  (5.25),  not  only  that  identical  boundary  values 
should  be  used  at  the  provisional  and  the  final  steps,  but  also  that  the 
mesh  size  must  be  sufficiently  refined  so  that  Re^x  < 2.  Sample  calcula- 
tions for  the  steady  state  solutions  of  the  linearized  Burgers'  equation 
(3.13)  verified  the  abrupt  change  of  the  behavior  from  a smooth  to  a 
violently  oscillatory  limiting  solution  when  Re^x  increases  beyond  the 
critical  value  of  2. 
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For  linear  problems  with  variable  coefficients,  the  various  funda- 
mental solutions  of  the  difference  equations  cannot  be  displayed.  It 
is  nevertheless  expected  that  the  spurious  solutions  will  be  suppressed 
if  the  same  operators  LjfT)  and  L2(T)  and  the  same  boundary  values  sre 
used  for  the  successive  iterative  steps  in  every  time  interval.  As  to 
the  proper  fundamental  solutions  of  [Lj(T)  ♦ l^OOjlK  ■ 0,  it  is  known 
that  one  of  them  must  be  unity  because  of  the  consistency  requirement. 

The  other  will  become  oscillatory  for  too  large  a Whether  the 

critical  value  of  Re^x  will  be  2 or  how  it  may  vary  with  x is  uncertain. 
For  nonlinear  problems  with  sufficiently  smooth  solutions,  the  completo 
suppression  of  the  spurious  fundamental  solutions  in  the  first  variation 
of  the  nonlinear  difference  operator  at  each  time  step,  may  be  expected. 
This  is  because  the  spurious  fundamental  solutions  contained  in  the  com- 
puted results  of  the  nonlinear  equations  will  have  been  reduced  to  higher 
order  small  quantities  in  At  by  the  stratagem  described  above.  Such- 
highor  order  small  quantities  in  At  will  be  of  little  significance  in  the 
steady  state  limit.  Ihus,  the  outstanding  problem  for  eliminating  shock- 
induced  oscillations  is  to  satisfy  the  requirement  of  using  some  suffi- 
ciently small  mesh  size  Ax  corresponding  to  the  restriction  of  Re^x  < 2 
for  the  linearized  Burgers'  equation.  It  is  anticipated  that,  for  non- 
linear problems,  there  may  not  be  such  a sharp  value  of  the  critical 
Re....  The  transition  from  smooth  to  oscillatory  steady  state  solution 

AX 

may  be  gradual  over  some  range  of  values  of  ReAx.  This  has  also  been 
verified  with  actual  computation.  It  is  supposed  that  the  following 
heuristic  model  will  give  a general  idea  where  this  critical  range  of 
Re^x  may  be' 
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When  • second  order  accurate  conservative  difference  algorithms  of 
the  class  (5.18)  is  used  for  the  integration  of  the  Navier-Stokes  equa- 
tions and  when  the  stratagem, just  described,  is  followed  in  the  treatment 
of  the  boundary  conditions,  the  shock  wave,  if  present  in  the  confutation, 
is  not  regarded  as  a discontinuity  but  a "smooth"  region  with  large  gra- 
dient, spread  out  by  the  pseudo-viscosity.  This  shock  transition  region 
usually  spreads  out  over  two  or  more  mesh  points  to  connect  the  smooth, 
asynptotically  uniform  flow  fields  both  up  and  downstream  of  the  shock 
region.  The  transition  profile  as  calculated  is  not  intended  to  be  accu- 
rate. Its  primary  function  is  to  accomplish  a smooth  connection,  and 
hopefully,  without  inducing  oscillations  propagating  into  the  smooth 
flaw  field  in  its  neighborhood.  Thus,  transition  profile  joining  a scalar 
function  u with  asymptotic  values  ua  > t a in  the  up  and  the  downstream 
region  respectively,  might  as  well  be  confuted  approximately  based  on 
the  nonlinear  Burgers'  equation  as  a model  for  the  local  flow  field.. 

This  means  that  the  local  profile  might  he  approximated  by  the  steady  solu- 
tion, Equation  (5.10)  with  x ■ 0 and  u ■ 0 located  at  the  point  of  maximum 
slope  in  the  transition  profile  that  has  been  actually  computed  with  the 

— dll 

full  Navier-Stokes  equations.  Thus,  the  computed  maximum  value  of 
properly  nondimensionalized  in  the  transition  region  will  define  the 
effective  Reynolds  number  of  the  transition  region. 

■ Re/2  (5.30) 

max  confuted 

In  this  manner,  the  poorly  defined  thickness  of  the  transition  region 
is  avoided.  The  parameter  a can  be  taken  as  imity  when  the  reference 
velocity  is  adopted  as  the  change  of  the  velocity  (or  the  particular 


m 


- 99  - 


scalar  quantity  in  dimensionless  form)  from  the  point  of  maximum  grad- 
ient to  the  asyiq>totic  value.  If  the  computed  transition  profile  is 
approximately  symmetric  with  respect  to  the  inflexion  point,  this  refer- 
ence velocity  will  be  half  the  juiap  across  the  shock. 

Let  the  asy^>totic  values  of  u across  the  shock  transition  region 


be  Uj  and  then,  assuming  > U2, 


/au\  . 

\ x / max  \ / max 


Ul.U2 


Hew  it  is  pressuwed  that  the  critical  value  of  this  Re^x  is  essentially 
the  same  as  if  the  computation  were  done  with  the  same  algorithm  but 
based  on  the  Burgers'  equation  so  that  oscillation-t.ee  confuted  results 
of  the  transition  region  can  be  effected  with  ReAx  < 2.  when  expressed  in 
terms  of  quantities,  directly  available  in  the  computation  as  an  a 
posteriori  criterion,  according  to  (5.31),  this  condition  becomes 


(AU)nn*  , 1 

r r ' 

i.e.,  "the  maximum  permissible  change  of  U per  mesh  (AU)^^.  is  one  half 
the  jus|>  |U1-U2I  across  the  discontinuity  so  as  to  avoid  shock-induced 
large  oscillations  in  the  computed  results." 

This  statement  implies  that  we  cannot  expect  to  obtain  an  oscilla- 
tion-free shock  front  that  contains  less  than  two  meshes  from  the  compu- 
tational solution  following  the  present  stratagem.  Moreover,  within  the 
linearized  framework,  the  criterion  (5.32)  should  be  equally  applicable 
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to  any  physical  scalar  variable,  sustaining  a "jiaqp"  across  some  large 
gradient  region  not  necessarily  a discontinuous  front,  although  the  Re^x 
was  defined  in  terms  of  flow  velocity  and  viscosity.  It  is  independent 
of  viscosity  explicitly. 

The  criterion  (S.32)  stands,  however,  only  as  an  a posteriori  cri- 
terion for  achieving  an  oscillation-frce  shock  solution.  This  is  because 

(AU)Bax  « f Ax  becomes  known  only  after  the  completion  of  the  compu- 

V /max 

tations;  by  then,  there  is  no  need  of  the  criterion  to  find  out  if  the 

computed  solution  is  oscillation-free!  Such  an  a posteriori  criterion 

can,  however,  be  of  some  help  in  practice,  since  (AU)  can  be  estimated 

max 

long  before  the  computed  solutions  readies  a satisfactory  "steady  state". 
Oscillations  will  be  present  in  the  "transient  states"  of  the  confutation 
whether  or  not  the  steady  state  limit  will  contain  shock-induced  oscilla- 
tions. If  the  criterion  should  be  satisfied  at  some  transient  stage,  we 
■ay  expect  an  oscillation -free  steady  state  solution  with  further  temporal 
steps.  Otherwise,  smaller  mesh  sizes  may  be  needed. 

It  is  more  convenient  if  this  criterion  can  be  put  into  some  a priori 
for*,  less  precise  as  it  must  be.  Note  that  the  magnitude  depends 

on  die  shock  strength,  the  shock  orientation  relative  to  the  coordinate 
axes  (in  a multi-dimensional  problem)  and  the  coordinate  direction  under 
consideration.  If  it  is  possible  to  estimate  this  magnitude  lUj-UjI,  then 

Re^-  |Ui-IJ2|Ax/2v  c 2 

may  be  used  directly  as  an  a priori  limit.  This  Reynolds  number  Re^x  must 
not  be  confused  with  the  Re^x  M based  on  the  uniform  supersonic  flow 
volocity  far  upstream  of  the  flow  field,  i.e.,  Re^x  m - 


U Ax/v.  In 


- 101  - 

tens  of  this  Re.„  . the  criterion  becomes 
AX,® 

U«A*  U.  4 

^Ax,®  “ v < l-Uj/Uj  (5.33) 

which  can  be  useful  epriori  if  there  is  sons  idea  as  to  the  shock  strength 

Uj/Uj  and  the  ratio  1 1,,/Uj  of  the  reference  velocity  far  upstream  to 

the  velocity  Uj  into  which  the  shock  wave  is  propagating.  For  complicated 

flow  problems,  however,  such  quantities  are  usually  among  the  unknowns. 

Thus,  the  limit  on  Re.  w given  by  (5.33)  will  have  to  be  based  on  some 

' , 

rough  estimate  or  on  the  "transient  states"  of  the  computed  solution. 

The  previous  heuristic  development  is  equally  applicable  to  any 
flow  region  containing  large  gradients  other  than  the  shock  front.  In 
particular,  oscillations  originating  from  boundaries  of.  the  field  of 
computation  can  be  likewise  alleviated.  It  is  to  be  emphasized,  however, 
that  if  the  oscillatory  extraneous  fundamental  solutions  like  5,  and 

O 4 

were  not  suppressed  by  the  stratagem  described  above,  these  extraneous 
oscillatory  solutions  will  propagate  into  the  neighboring  smooth  flow  fields 
even  if  the  mesh  size  is  much  reduced  below  what  is  required  by  (5.33),  at 
least  one  of  them  will  be  amplifying  away  from  the  boundaries  of  the  trans- 
ition region  while  propagating  into  the  neighboring  smooth  regions  on 
either  side.  On  the  other  hand,  if  much  too  coarse  a mesh  size  is  used  in 
the  computation,  large  amplitude  oscillations  will  result  since  one  of  the 
proper  fundamental  solutions  of  the  difference  equation  fails  to  be  a valid 
approximation  to  that  of  the  differential  problem  despite. the  fact  that  the 
stratagem  described  above  is  followed.  To  produce  an  os  dilation -free  com- 
putational solution  of  a flow  problem  involving  shock  waves,  it  is  recom?- 
mended  not  only  that  some  form  of  the  two  step  algorithm  (5.18)  be  used 
with  identical  botndary  values  applied  to  both  iterative  steps  during  a 
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tints  interval,  but  also  that  the  mesh  size  Ax  is  kept  sufficiently  small 
compared  with  the  condition  (5.33).  This  recommendation  is  based  on  the 
results  of  analysis  of  a simple  model  linear  equation  for  the  numerical 
solution  of  the  much  more  complicated  and  nonlinear  gas  dynamic  equations. 
It  is  recommended  in  the  some  spirit  that  the  local,  linear  stability 
analysis  of  v. Neumann  will  help  in  achieving  computational  stability. 

The  practical  merit  of  this  recommendation  is  yet  to  be  examined  in 
greater  detail  by  the  computational  community. 

The  previous  development  has  guided  the  author  quite  successfully 
in  his  earlier  attempts  of  integrations  of  the  N&vier-Stokes  equations 
for  sore  complicated  flow  problems,  such  as  the  near  wake  flow  behind  a 
flat  base  with  a sharp  comer  in  a supersonic  floi'Ud  the  hypersonic 
flow  over  the  sharp  leading  edge  of  a highly  cooled  flat  plate^lhe 
flow  situations  encountered  in  these  examples  are  just  too  complicated 
to  provide  any  meaningful  quantitative  tests  of  the  validity  of  this 
criterion  and  the  accuracy  of  the  confuted  results.  In  the  following 
a simple  case  will  be  described  which  may  serve  to  support  and  to  illus- 
trate the  usefulness  of  the  stratagem  and  the  simple  criterion  despite 
the  heuristic  content  of  its  application  to  the  actual  integration  of  the 
Navier-Stokes  equations. 

The  Gheng-Allen  two  step  algorithm  as  a member  of  the  class  (5.18) 
is  used  to  integrate  the  complete  Navier-Stokes  equations  for  the  propa- 
gation of  a planar  shock  wave  into  a uniform  supersonic  flow  at  Mach  No. 

2 with  the  shock  front  inclined  at  an  angle  8 « 41.84°  to  the  uniform  in- 
gas  density  Pj,  velocity  Uj, energy  «j,  and  pressure  pj  are  taken 
to  be  unity  in  dimensionless  form.  The  theoretical  values  of  these  vari- 
ables downstream  of  the  shock,  according  to  the  Hugoniot  relations  agree 


with  the  values  computed  at  Re^x  w • 10  to  better  than  0.1%.  The  criti- 
cal Reynolds  nuaber  per  mesh  is  (ReAx  ^ )c  ■ 4/(1-0.837)  « 24.5.  No 
oscillations  are  found.  The  shock  front  is  sharp  and  straight.  It  is 
verified  that  the  a posteriori  criteria  (5.32)  are  satisfied  for  the 
density  p the  x-velocity  component  u,  the  y-velocity  component  v,  the 
energy  e and  the  pressure  p across  the  shock. 

When  the  computation  is  repeated  at  Re^x  m • 50  exceeding  the  criti- 
cal value  (ReAx-)c  * 24.5  for  the  same  flow  configuration,  sifcstantial 
oscillations  are  present  immediately  downstream  of  the  shock.  The  a 
posteriori  criteria  (5.32)  for  all  the  physical  variables  are  found  vio- 
lated. The  peak  amplitude  of  the  oscillation  is  about  10%  but  such  oscil 
1 atioas  are  essentially  damped  out  a few  meshes  downstream  of  the  shock. 
The  downstream  asym>totic  values  are  reached  well  within  the  field  of 
computation.  The  downstream  asymptotic  values  obtained  from  the  computa- 
tion at  ReAx  m • 50  are  correct  to  within  0.3%  of  the  Hugoniot  values. 

The  smooth  incident  shock  computed  at  ReAx  ^ ■ 10  was  then  allowed 
to  be  reflected  from  an  inviscid  wall.  For  the  reflected  shock,  the  crit 
ical  Reynolds  number  (Re^  )£  » 4/(0.837  - 0.646)  * 21,  which  exceeds 
the  ReAx  m ■ 10  used  in  the  computation.  A smoc  •t»,  straight  reflected 
shock  is  obtained.  All  the  computed  downstream.  iptotic  values  agree 
with  the  theoretical  values  to  bette.  than  0.1%  and  there  are  no  oscilla- 
tions. 

Computations  at  intermediate  values  of  Re^x  indicate  that  oscilla- 
tions begin  to  appear  with  ReAx  exceeding  10  to  15(  increase  most  rapidly 
around  the  critical  value  of  20  - 30,  and  keep  increasing  slowly  with 
larger  ReAx.  This  gradual  rather  than  an  abrupt  change  of  behavior  with 
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probably  what  should  be  expected  in  a nonlinear  system.  It  is 
encouraging  that  the  sicple  criterion  obtained  from  an  elesentaxy  linear 
analysis  of  a simple  model  may  prove  to  be  useful  in  complicated  flow 
problems  to  be  encountered  in  practico. 


VI.  CURRENT  STATUS  AND  FUTURE  PROSPECT 


The  various  problems  associated  with  the  numerical  integration  of 
Navier-Stokes  equations  have  been  reviewed  in  the  previous  chapters  as 

to  the  mathematical  origin  of  the  problems  and  of  the  basis  of  various 

/ 

techniques  in  dealing  with  them.  This  approach  was  chosen  in  prefer- 
ence to  a review  in  the  form  of  a glossary  of  various  solutions  in  the 
literature  so  as  to  provide  a frame  of  reference  how  such  solutions 
may  be  studied  and  how  each  specific  problem  should  be  approached. 

In  the  days  of  the  mechanical  desk  calculators  or  the  card  pro- 
grammed calculators  (CPC),  the  numerical  integration  of  the  hydrodynamic 
equations  was  attempted.  The  primary  concern  then  was  the  limitation 
of  the  computational  speed  offered  by  these  mechines.  While  the  question 
of  cos^wtational  stability  was  known  to  mathematicians^,*  it  is  not  of 
much  concern  to  the  practitioners.  The  dawn  of  the  high  speed  electronic 
computers  in  the  mid-1940's  changed  all  that.  The  ability  to  compute  fast 
showed  how  often  an  apparently  straight forward  computation  will  lead  to 
unbounded  meaningless  results.  This  problem  is  the  first  and  the  most 
pressing  one  presented  by  the  high  speed  computers.  If  the  stability  ques- 
tion of  the  computation  is  not  successfully  resolved,  no  results  of  any 
kind  could  be  obtained.  Since  the  mid  1940's,  this  stability  question  has 
been  studied  very  extensively,  both  mathematically  and  empirically.  As 
was  described  in  Chapter  III,  much  has  been  learned  and  understood  since 
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then.  It  is  still  true  that  when  complicated  sets  of  partial  differen- 
tial equations  such  as  those  of  gas  dynamics  are  to  be  integrated,  compu- 
tational  stability  remains  a formidable  problem.  As  in  the  older  days, 
so  much  work  is  still  needed  to  render  a stable  computation  that  one 
often  hesitates  to  ask  any  further  questions  about  any  reasonably  looking 
computed  solution.  But  for  those  interested  in  cooputational  methods  for 
some  practical  purposes,  computational  stability  is 'not  synonymous  with 

i 

the  major  problem  of  the  computational  solution  of  a partial  differential 
equation  system.  It  is  only  a first  step  in  achieving  a solution  of  value 
in  practice. 

With  the  help  of  suitable  model  studies  and  appropriate  choices  of 
difference  algorithms,  computational  stability  can  generally  be  obtained 

and  tested  in  actual  machine  computation.  Now  is  the  time  to  be  concerned 
with  obtaining  not  only  some  qualitatively  correct  solutions  but  also 
quantitatively  correct  answers  with  some  estimate  of  the  error  bounds  of 
Idle  computed  solution.  In  applications,  the  primary  purpose  of  a computed 
solution  is  to  seek  some  reasonably  accurate  quantitative  estimate  of  the 
flow  field.  The  accuracy  requirements  for  different  applications  may  be 
quite  different.  Whether  a solution  is  sufficiently  accurate  for  a spe- 
cific application  can  only  be  judged  under  some  overall  view.  But  such  a 
judgement  can  be  made  only  when  the  computed  solution  is  accompanied  by 
an  error  bound  if  not  a strict  error  estimate.  The  error  bounds  of  a 
computed  solution  is  no  less  important  than  the  error  bars  of  a set  of 
experimental  data  if  such  computed  or  experimental  sets  of  data  are  to 
be  practically  useful.  With  this  in  view,  the  preliminary  developments 
on  computational  accuracy  given  in  Chapters  IV  and  V are  very  important 
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in  practice.  Most  of  the  solutions  available  in  published  literature 
cannot  be  examined  with  regard  to  the  question  of  accuracy  because  they 
were  probably  obtained  primarily  to  demonstrate  qualitatively  what  can 
be  doen  rather  than  to  solve  specific  problems  in  application.  A few 
examples  will  be  described  below  with  some  comment. 

6.1  Hydrodynamics 

The  flow  of  an  incompressible  viscous  fluid  in  two  space  dimension 
probably  represents  the  simplest  form  of  the  Navier-Stokes  equations. 

It  is  most  often  treated  in  the  stream-function-vorticity  form.  The 
mass  continuity  equation  in  two  space  dimension  (x,y) 


3u  ^ 3v 
37  5y 


0 


(6.1) 


can  be  satisfied  by  a scalar  stream  function  V defined  with 

M . 3Y 

u " 5y  and  v “ - 57  (6.2) 

while  the  vorticity  component  o>  normal  to  the  x-y  surface  is 


V2Y 


92V  d2V 

5*r  ♦ 
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(6.3) 


The  curl  of  the  momentum  equation  reduces  to  the  vorticity  transport 
equation 


3w  St  3w  9f  3b) 
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(6.4) 


The  divergence  of  the  momentum  equation  gives  the  V*p  in  terms  of  ¥ and 
u.  Thus  the  static  pressure  in  the  flow  field  can  be  solved  independently 
after  the  stream  function  f and  vorticity  ca  have  been  determined.  Thus 
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the  solution  for  a hydrodynamic  problem  will  be  posed  as  the  solution  for 
two  elliptic  problems  for  4*  and  u simultaneously  represented  by  Equations 
(6.3)  and  (6.4)  subject  to  Dirichelet  and/or  Neumann  boundary  conditons 
on  a closed  boundary  for  ¥ and  tu.  The  physical  bouidary  condition  depends 
on  the  problem. 

A simplest  case  is  the  decay  of  a vortex  in  a closed  rectangular  box 
in  which  case  u • v ■ 0 on  the  boundary  which  stay  be  taken  as  x « 0,  y * 0, 
x » 1,  y - 1,  respectively.  This  set  of  physical  boundary  conditions  has 
to  be  translated  into  boundary  conditions  of  ¥ and  u.  By  definition  V - 0 
may  be  assigned  on  the  boundary.  This  serves  to  determine  ¥(x,y)  completely 
from  (6.3)  when  w(x,y)  is  given  over  the  field.  The  remaining  physical 
boundary  conditions  are 

2_3u*0  on  x * 0 or  x * 1 

BY 

- ■ v » 0 on  y » 0 or  y - 1 (6.5) 

The  practical  question  arises  how  (6.5)  may  be  expressed  as  the  boundary 
conditions  of  at  in  the  solution  of  Equation  (6.4).  In  practice,  this 
question  is  by-passed  by  solving  Equation  (6.3)  first  for  the  advanced 
values  of  ¥(x,y)  and  the  boundary  values  of  0)  are  estimated  either  from 
the  initial  data  or  the  results  of  the  most  recently  available  advanced 
values  of  f near  the  boundary.  This  can  be  done  either  with  or  without 
the  conditions  (65)  taken  into  consideration.  In  principle  the  boundary 
conditions  (6.5)  should  at  least  be  checked  a posteriori.  There  is 
clearly  an  error  tg  on  the  boundary  values  of  u>  of  the  order  of  At,  Ax, 
and/or  Ay  depending  on  the  formal  order  of  accuracy  how  this  boundary  con- 
dition is  handled. 
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Now  if  Equation  (6.4)  is  integrated  over  the  volume  (i.e.,  x ■ 0 to 
1 and  y • 0 to  1)  and  over  the  time  period  t « 0 to  t of  the  integration, 
the  total  decay  of  the  vorticity  is 


f lu»0(t  = 0)  - u(t))<V 

•'v 

' v/ dt/  72oj  dV 


i.e.,  proportional  to  the  total  outflux  of  the  gradient  of  vorticity  on 
the  boimdary.  (The  three  dimensional  analog  is  obvious).  Thus  the  non- 
random cumulative  error  on  the  total  decay  of  the  vorticity  in  the  box 
will  be  of  the  order  of  NJXg  where  N is  the  number  of  time  steps  inte- 
grated and  J is  the  number  of  spatial  meshes  in  a linear  dimension.  The 
$ • 0)  on  the  boundary  is  assumed  to  be  of  the  same  order  as  the  error  in 
the  boundary  vorticity  xB  itself.  The  use  of  the  integral  formula  implies 
that  the  accumulation  of  the  truncation  errors  over  all  the  interior  points 
have  been  neglected  in  the  difference  formulation.  Even  so,  the  total 
decay  of  the  vorticity  at  later  times  depend  very  importantly  on  how 
accurately  the  boundary  vorticity  was  formulated  in  the  computation  end 
on  whether  and  hew  the  errors  associated  with  such  a formulation  will  accum- 
ulate in  space  (along  the  boimdary  ) and  in  time.  The  question  is  more 
than  the  local  truncation  error  of  the  difference  formulation  of  the  vor- 
ticity boundary  condition  since  the  correct  physical  boundary  condition 
Equation  (6.5),  which  represents  some  integrated  condition  on  the  vortc- 
ity  field  rather  than  the  local  values  of  the  vorticity  was  ignored. 

The  me  of  the  stream  function -vorticity  as  the  dependent  variables 
ia  the  fundamental  reason  of  the  difficulty  in  implementing  the  boundary 
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condition.  It  also  causes  considerable  conq>li cations  in  rendering  a 
conservative  formulation  to  prevent  the  accumulation  of  the  truncation 
errors  over  the  interior  points.  If  the  physical  variables  u and  v are 
used  as  the  dependent  variable  in  the  difference  foraulation,  the  diffi- 
culty with  the  boundary  cc.  Jition  would  be  eliminated  for  the  above  ex- 
ample and  the  conservation  of  the  difference  formulation  can  be  readily 
implemented.  Ihe  advantage  of  the  vorti city -stream  function  formulation 
in  reducing  the  nuober  of  partial  differential  equations  may  be  more  than 
compensated  by  the  difficulty  it  brings  for  this  problem. 

For  hydrodynamic  problems  with  inflow  and  outflow  boundary  in  the 
field  of  corap utati on,  the  boundary  treatment  in  the  difference  formula- 
tion faces  a difficulty  of  different  nature.  This  is  because  the  physi- 
cal boundary  conditions  are  prescribed  very  far  up  and  downstream  of  the 
field  of  computation.  The  vorti city-stream  function  formulation  does  not 
aggravate  the  situation  much  further  and  is  therefore  often  preferred  for 
the  numerical  integration  of  the  hydrodynamic  equations.  The  Poisson  type 
equations  can  be  efficiently  solved  in  different  ways.  There  are  many 
such  solutions  in  the  literature.  Most  of  such  results  can  not  be  analyzed 
for  an  error  estimate  primarily  because  of  the  non- conservative  form  of 
the  difference  formulation  that  permits  the  accumulation  of  the  local  trun- 
cation errors.  Experimental  data  are  generally  not  available  to  provide 
a quantitative  estimate  of  the  error  of  the  computed  results.  All  sudi 
computations  serve  to  demonstrate  the  feasibility  of  confuting  some  "reason- 
able" approximate  solutions  but  are  of  little  quantitative  value.  A numer- 
ical study  of  the  steady  flow  of  a uniform  stream  over  a sphere  was  there- 
fore undertaken  by  Rimon  and  Qieng. 


The  flow  field  of  a uniform  stream  over  a sphere  is  conveniently 
described  by  using  the  spherical  polar  coordinates.  To  extend  the  outer 
boundary  of  the  field  of  coogmtation  to  os  far  downstream  as  possible  to 
facilitate  the  implementation  of  the  boundary  conditions,  i ■ In  r is 
used  in  place  of  the  physical  radius  r.  Three  different  sets  of  numeri- 
cal integration  have  been  made  by  different  authors  at  common  Reynolds 
numbers  of  40  and  100.  There  is  also  a set  of  experimental  data  by  Taneda 
of  some  characteristic  quantities  of  the  recirculatory  wake  flow  field 
at  these  end  other  Reynolds 'numbers.  Such  measured  values  of  wake  length, 
locations  of  the  separation  point  and  the  vorticity  centers  provide  com- 
parisons of  the  detailed  flow  field  in  the  most  sensitive  region,  in  addi- 
tion to  the  overall  drag  coefficient  acting  cn  the  sphere. 

[20]  [21] 

Jenson  and  Hamielec,  et  al.  used  similar  difference  relaxation  pro- 
cedures and  used  the  same  downstream  boundary  conditions  approximating 
imiform  out  flow.  Both  cases  were  carefully  executed  and  examined  numer- 
ically, very  carefully,  and  made  sure  that  the  steady  state  results  they 
obtained  are  essentially  independent  of  further  reduction  of  mesh  spacing 
from  mesh  sizes  A8  ■ 6*  and  z * 1/20.  They  obtained  the  same  drag  coef- 
ficient CD  in  agreement  with  what  is  expected  from  the  experimentally 
well-established  standard  drag  curve.  However,  the  details  of  the  two 
solutions  were  much  different.  For  example,  the  vorticity  on  the  wake 
side  of  sphere  surface  differ  by  a factor  of  2 to  3 for  the  case  with 
Rep  » 40.  The  streamline  patterns  in  the  recirculatory  wake  are  visibly 
different  although  qualitatively  similar.  It  is  supposed  that  such 
differences  in  the  results  largely  demonstrate  the  cumulative  effects  of 
the  truncation  errors  due  to  the  non -conservative  nature  of  the  difference 
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algorithms  equivalent  to  their  relaxation  procedures.  Jenson's  results 

depart  considerably  further  from  Taneda's  wake  data  than  the  results  of 

/ 

Hamielec,  at  al.  at  ReQ  ■ 40.  Hamielec,  et  al.  also  calculated  the  case 

ROq  » 100.  They  found  it  necessary  to  refine  the  mesh  to  A0  « 3s  and 

Az  « 1/40  to  secure  a reasonable  steady  state  and  to  introduce  some  fine 

adjustments  in  order  to  reproduce  the  experimental  value  of  the  drag 

coefficient  C_  at  ReQ  ■ 100. 

[22] 

Risen  and  Cheng  succeeded  in  developing  a conservative  difference 

fora  that  is  still  reasonably  simple  despite  the  contracted  curvilinear 

coordinates  and  stream  function  vorticity  formulation.  The  same  mesh 

size  A0  - 6°  and  Az  - 1/20, as  was  used  by  the  previous  authors,  was  used. 

The  conservative  nature  of  the  difference  formulation  permit  an  estimate 

of  the  upper  bound  of  the  cumulated  truncation  error  by  Equation  (5.14). 

The  Re^x  should  be  replaced  by  Re^z  for  this  calculation  in  terns  of 

A9  and  Az.  The  magnitudes  of  Re^z  for  the  two  cases  with  Re^  * 40  and 

100  can  be  estimated  from  the  computed  solution  based  upon  the  velocity 

gradient  in  the  region  near  the  isolated  rear  stagnation  point  in  the 

wake  (the  reference  velocity  of  this  mesh  Reynolds  number  is  based  on 

the  velocity  difference  across  the  large  gradient  region).  They  are  less 

than  j end  1 respectively.  Accordingly,  the  absolute  uq>per  bound  of  the 

-2  2 

cumulated  truncation  errors  are  3 x 10  x Re^z  % I end  3%  respectively. 
Ihe  extrapolation  condition  at  the  downstream  boundary  gives  rise  to  the 
largest  contribution  to  the  boundary  error.  Both  ^ 0 and  u > 0 on  the 

out  flow  boundary  may  commit  a fractional  error  as  much  as  100%.  (It  is 
not  ejected  to  error  in  sign).  The  absolute  upper  bound  of  the 

i 

boundary  errors  may  then  be  estimated  with  Equation  (5.16)  where  * 1 
and  Re  is  based  on  the  maximum  flow  velocity  in  the  wake  region  and  the 
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length  from  the  rear  stagnation  point  to  the  out  flow  boundary  in  the 
z-0  plane  of  computation.  This  is  more  than  two  unit  lengths.  The  ef- 
fective Reynolds  numbers  are  then  80  and  200  respectively.  Accordingly, 
the  bounds  of  the  boundary  errors  are  estimated  as  2e^/Re  of  2.5%  and  1% 
respectively.  By  adding  the  estimates  of  the  absolute  upper  bounds  of 
the  truncation  errors  and  the  boundary  errors  for  each  case,  the  overall 
estimate  of  the  absolute  error  bounds  are  about  3.5%  and  4%  for  the 
cases  ReD  * 40  and  100  respectively.  This  is  a quite  satisfactory  engin- 

t 

eering  accuracy.  Thus,  it  is  expected  and  verified  that  the  computed 
results  of  Rimon  and  Cheng  will  agree  with  Taneda's  wake  data  much  better 
than  the  results  of  Jenson  end  Hamielec,  et  al.  The  computed  vorticity 
field  in  the  near  wake  region  of  Rimon  and  Cheng  and  that  from  Hamielec, 
et  al.  differ  by  a factor  of  2 or  more  in  the  case  with  Re^  = 100  while 
they  differ  much  less  for  the  case  Re^  = 40.  This  again  demonstrates 
the  significance  of  the  accumulation  of  the  local  truncation  errors. 

The  computational  effort  of  the  solution  of  this  problem  following 
the  formulation  of  Rimon  and  Cheng  was  not  excessive  at  the  time  and  is 
rather  small  in  terms  of  the  present  day  computing  machines.  61  x 31  * 

1691  mesh  points  were  used.  Steady  state  solution  was  obtained  in  about 
an  hour  computation  per  case  in  the  IBM  7094,  with  the  potential  flow 
as  the  initial  data.  In  terms  of  CDC  6600  machine  time,  it  would  take 
less  than  10  minutes.  The  computational  time  can  be  appreciably  reduced 
if  a more  reasonable  approximation  than  the  inviscid  flow  field  should  be 
used  as  the  initial  data.  It  is  therefore  believed  that  with  conscientious 
effort  in  constructing  the  difference  formulation,  useful  quantitative 
results  can  be  obtained  from  the  numerical  solution  of  the  Navier-Stokes 
equations. 
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The  extension  of  such  calculations  to  steady  flows  in  three  dimensional 
space  and  to  higher  flow  Reynolds  numbers  will  need  not  only  substantially 
■ore  computer  tine  but  also  some  analysis  to  gain  understanding  of  certain 
intricacies  in  the  3-D  problems  especially  concerning  the  computations  in 
the  vicinity  of  "separation  lines".  Kith  the  greatly  increased  capability 
of  the  high  speed  computers  in  the  foreseeable  future,  it  is  reasonable  to 
expect  good  quantitative  results  for  these  steady  3-D  problems  of  practical 
interest. 

/ 

The  time  dependent  hydrodynamic  problems  in  three  space  dimensions  may 
be  considerably  more  difficult  and  demanding.  This  is  especially  true  if 
the  hydrodynamic  turbulence  is  the  subject  of  investigation.  The  high  fre- 
quency components  of  the  turbulent  fluctuations  could  doubtfully  be  treated 
with  a reasonable  accuracy  despite  the  giant  stride  in  the  capability  of 
the  computing  machines  foreseen  in  the  future.  It  appears  that  some  phe- 
nomenological theory,  at  least  for  the  high  frequency  components,  will  be 
needed  while  the  low  frequency  components  may  be  satisfactorily  handled 
by  computational  methods.  This  statement  is  meant  to  apply  whether  it  is 
to  be  integrated  in  the  physical  space  with  the  physical  variables  or  it 
is  to  be  treated  in  the  Fourier  space  for  the  Fourier  components  of  the 
physical  variables.  Much  work  is  needed  in  any  case f2**25, 26 ^ 

6.2  Supersonic  Gas  Dynamics 

The  gas  dynamic  equations  system  is  basically  the  same  as  the  hydro- 
dynamic  equation1-  . xcept  for  the  variati  is  of  gas  density  and  the  dif- 
fusivities  and  for  the  addition  of  the. equation  of  energy  balance  (1.3). 

The  outstanding  feature  of  the  supersonic  flow  field  is  the  presence  of 
shock  waves  either  generated  from  within  or  incident  on  the  flow  field 
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from  without.  Most  of  the  practical  problems  that  call  for  the  numerical 
treatment  of  the  Navier  Stokes  equations  involve  the  generation  of  shock 
wave  froD  the  interaction  between  the  inviscid  and  the  viscous  streams. 

The  computation  of  a shock  wave  with  unknown  strength  and  location  presents 
considerable  difficulties  as  was  discussed  in  the  previous  chapter.  The 
shock-induced  oscillations  in  its  neighboring  flow  field  is  detrimental 
to  the  appearance  of  the  computational  solution.  Such  solutions  were 
often  presented  after  some  artificial  averaging  or  filtering  procedure  and 
can,  therefore,  be  of  qualitative  values  only.  Those  solutions  relatively 
free  frea  this  criticism  indeed  ewe  their  success  in  avoiding  the  serious 
consequences  of  a shock  standing  in  an  important  part  of  the  flow  field* 

By  carefully  selecting  the  field  of  computation  for  the  problems  to  be 
investigated,  they  minimized  the  consequences  of  shock-induced  oscilla- 
tions . 

[16] 

Allen  and  Oieng  treated  the  near  wake  flow  inbedded  in  a supersonic 
stream  turning  over  a sharp  shoulder  with  a "recompression  shock"  generated 
from  the  turning  of  the  supersonic  stream  caused  by  the  closing  of  the 
recirculatory  wake.  In  the  steady  state  solution  of  this  problem,  the 
small  oscillations  caused  by  the  recompression  shock  distorts  appreciably 
the  computed  results  only  in  the  far  downstream  portionof  the  rejoined 
wake  flow  field  near  the  downstream  boundary.  Although  the  oscillations 
of  the  flow  properties  in  the  flow  field  is  equivalent  to  those  induced 
by  an  oscillation  of  the  shock  front  of  only  l/4Ax,  it  .remains  as  or.e  of 
the  two  largest  sources  of  computational  errors.  It  is  conjectured  that 
the  likely  source  of  the  small  oscillation  is  the  extraneous  inaccurate 
difference  treatment  where  the  shock  emerges  from  the  downstream  boundary 
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of  the  field  of  computations.  The  conservative  difference  form  of 

the  class  (5.18)  was  used  and  the  criterion  (5.33)  is  satisfied  although 

without  a substantial  margin.  Unfortunately,  comparable  experimental 

data  are  not  available  and  the  extension  of  this  calculation  to  the  range 

3 4 

of  practical  Reynolds  numbers  of  10  - 10  and  for  the  somewhat  more  com- 

plicated geometrical  configuration  was  beyond  the  means  (cost  of  computa- 
tions) available. 

[27] 

Ross  and  Cheng  studied  the  question  that  if  the  number  of  mesh  points 
is  limited  to  2100  with  "optimal"  ratio  of  Ax/Ay  and  with  some  non-essen- 
tial but  simplifying  modifications  of  the  boundary  treatment,  what  the  range 
of  the  Reynolds  numbers  and  Mach  numbers  will  be  when  the  computational 
solutions  with  the  previous  formulation  will  possess  an  absolute  upper,, 
bound  of  the  error  of  no  more  than  10%.  The  computational  effort  was 
limited  to  10-15  minutes  of  computing  time  on  the  IBM  360-91  equivalent 
roughly  to  20-30  minutes  on  the  CDC  6600  or  4 to  10  hours  on  the  IBM  7044, 
originally  used  by  Allen.  When  other  restrictions  purely  of  fluid  mechan- 
ical nature  are  superposed,  it  was  established  that  the  range  of  validity 
of  the  computational  formulation  can  be  extended  to  M = 4 and  Reynolds  num- 
bers of  'v*  1-2  * 10^  based  on  half  width  of  the  base.  To  extend  this  compu- 
tation to  the  practical  range  of  interest  would  require  substantial  re- 
' finenent  in  the  mesh  size  with  corresponding  increase  of  the  computational 
effort.  The  storage  requirement  on  the  computer  does  not  seem  restrictive. 
It  was  the  computer  time  and  cost  that  was  prohibitive. It  may  be  that  the 
absolute  upper  bound  of  10%  is  too  restrictive  since  the  maximum  fractional 
error  in  the  solution  is  likely  to  be  substantially  less  than  the  absolute 
upper  bound.  A substantial  decrease  of  the  estimate  of  the  computational 
effort  will  follow  a modest  reduction  of  the  accuracy  requirement  if  the 
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oothod  of  error  estimate  described  in  Chapter  V should  be  granted  in  the 
absence  of  any  direct  comparison  with  reliable  and  comparable  experimental 
data. 

Carter2 o^ooses  to  integrate  the  Navier-Stokes  equation  for  a steady 
supersonic  viscous  flow  over  a compression  ramp  or  comer  with  an  iabedded 
separated  region.  Ihe  compression  waves  will  eventially  coalesce  into  a 
shock  wave.  Carter  kept  the  upper  boundary  of  the  field  of  computation 
sufficiently  close  to  the  viscous  region  so  that  the  waves  generated  from 
the  viscous  layer  any  be  treated  as  isentropic  waves  without  serious  error 
and  utilized  the  sicple  wave  extrapolation  condition  on  the  upper  boundary. 
This  stratagem,  as  was  used  in  the  treatment  of  the  near  wake  problem, 
serves  to  eliminate  the  major  part  of  the  undesirable  wave  reflection  from 
the  iqpper  boundary.  By  restricting  the  field  of  confutation  to  such  a 
narrow  strip  and  using  a highly  refined  mesh  with  Brailovskaya's  difference 
algorithm,  a member  of  the  class  of  (5.18),  the  results  compare  favorably 
with  experimental  data  in  the  cotq> arable  Reynolds  and  Mach  number  ranges. 
The  difference  formulation  is  probably  not  quite  conservative  due  to  the 
use  of  the  "curved”  body  coordinates.  But  the  curvature  is  sufficiently 
smell  or  otherwise  localized  so  that  the  accumulation  of  the  truncation 
errors  may  not  be  excessive.  While  an  estimate  of  the  error  bounds  has 
not  yet  been  made,  the  evidence  seems  to  indicate  that  this  calculation 
may  have  come  very  close  to  generating  directly  some  useful  practical  re- 
sults. Admittedly,  the  computational  effort  in  this  calculation  seemed 
to  be  excessive  from  the  academic  point  of  view,  (two  or  more  hours  of  CDC 
6600  per  case) , it  does  not  appear  prohibitive  from  the  view  of  engineer- 
ing development  Moreover,  there  is  sii)stantial  room  for  improvements 
if  an  error  estimate  can  be  made.  The  4th  generation  computers  that  will 
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shortly  be  operational  further  promise  a substantial  increase  in  speed  of 
confutations  and  in  storage  capability.  This  may  render  the  computational 
effort  to  be  of  less  concern  in  solving  such  practical  problems. 

An  academic  program  has  been  devoted  to  developing  techniques  to  face 
the  difficulty  of  confuting  shock  waves  in  a complicated  flow  field.  Ihe 
results  reported  in  Section  5.4  demonstrated  some  progress  in  this  direc- 
tion. There  are  still  tremendous  difficulties  ahead  and,  as  yet,  not  re- 
solved whon  the  shock  wave  interacts  with  other  incident  waves  and  when 
the  criterion  of  (5.33)  becomes  much  too  restrictive.  Nevertheless,  even 
in  the  presont  unsatisfactory  state,  the  conf utationel  results  can  be 

useful  in  fluid  dynamics  research  to  supplement  experimental  and  other 

« 

efforts.  They  will  also  encounter  difficulties,  but  of  some  different 
nature.  The  following  treatment  of  the  hypersonic  leading  edge  problem 
may  illustrato  the  situation. 

Over  the  leading  edge  of  an  infinitely  thin  flat  plate,  placed  in  a 
hypersonic  or  supersonic  stream  at  zero  incidence,  a shoe*  wave  will 
develop  due  to  the  viscous  effects  in  the  vicinity  of  the  plate.  In  this 
region,  the  hypersonic  strong  interaction  theory,  based  on  the  boundary 
layer  type  arguments,  failed  to  provide  even  a qualitatively  adequate 
description  of  the  flow  field.  It  is  in  doubt  to  what  extent  the  flow 
situation  will  have  to  be  described  by  the  kinetic  theory  rather  than  by 
the  continuum  theory  when  appropriately  modified  for  tho  slip  effects. 

The  flow  field  was  thus  computed  by  Chen  and  Cheng.  A rather  strong 
oblique  shock  wave  develops  rapidly  from  the  leading  edge,  and  produces, 
in  the  downstream  gas,  a high  pressure  and  teiqperature , both  proportional 
to  Masin*0  where  6 is  the  local  inclination  of  the  shock  front  with  the 
incooing  uniform  stream.  It  is  very  clear  then  that  any  small  oscillations 
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in  the  shock  front  will  produce,  in  the  downstream  corresponding  oscilla- 
tions of  very  significant  magnitudes  with  the  upstream  flow  Mach  number 
of  *\*  20.  It  is,  therefore,  very  critical  to  have  the  shock  oscillations 
essentially  eliminated  from  the  computation.  The  conservative  two  stop 
algorithm  of  Cheng  and  Allen  was  again  used  with  a 40  x 30  mesh  in  the 
physical  space  x-y  with  y ■ 0 describing  the  plate  surface.  The  loading 
edge  shock  emerges  from  the  downstream  out  flow  foundary  just  below  the 
top  coiner.  No  oscillation  of  the  shock  wave  is  noted  except  in  the 
immediate  vicinity  of  the  plate  leading  edge  point.  The  oscillation  is 
fairly  large  but  dies  out  rapidly  within  about  S noshes  downstream  along 
the  plate  and  2-3  meshes  normal  to  the  plate.  Various  slip  conditions 
were  used  in  the  computation.  A minor  localized  oscillation  developed 
somewhere  downstream  about  2-4  Ay  for  no  obvious  reasons  and  is  far  away 
from  the  shock.  It  is  conjectured  to  have  originated  from  some  inapprop- 
riate treatment  of  the  boundary  conditions  on  the  plate.  This  localized 
oscillation  imposes  no  significant  error  on  the  solution.  The  downstream 
out  flow  boundary  is  again  treated  by  the  second  order  accurate  extrapola- 
tion but  along  the  shock  direction  where  the  shock  emerges,  along  the 
plate  cn  the  plate  surface  and  along  directions  linearly  interpolated  in 
between.  Despite  the  suppression  of  the  oscillations  of  the  emerging 
shock,  the  slow  decay  of  this  boundary  error  stands  as  the  largest  singlo 
contributor  to  the  solution  in  the  interior.  The  absolute  tqpper  bound 
is  <7%  according  to  Equation  (5.16)  as  evaluated  with  the  smooth  coqmted 
solution  a posteriori.  Tho  truncation  error  as  evaluated  with  Equation 
(S.14)  and  the  solution  in  the  non-os cillatozy  region  away  from  the  leading 
edge  point  is  less  than  7%.  With  both  the  routd-off  error  and  the  error 


due  to  the  steady  state  criterion  both  less  than  1%,  it  is  estimated  that 
the  absolute  upper  boind  of  the  error  in  the  cosputcd  solution,  away  from 
the  immediate  vicinity  of  the  leading  edge  and  the  out  flow  boundary,  is 
16%.  The  scatter  of  experimental  data  is  more  than  50%. 

Comparison  of  the  computed  results  with  a collection  of  experimental 
data,  not  precisely  corresponding  but  encompassing  the  case  computed,  shows 
sifcstantial  agreement  with  the  data  considered  to  be  most  reliable  in  the 
different  regions  and  for  different  quantities.  This  comparison  is  cer- 
tainly not  the  best  and  the  most  definitive,  but  is  probably  the  best  avail- 
able and  possibly  the  best  one  may  hepe  to  have  in  the  not  too  distant 
future.  This  is  in  view  of  the  difficulty  in  reproducing  the  experimental 
environments  and  of  the  cost  involved  in  such  experiments.  This  study 
leads  to  a few  physically  meaningful  conclusions  which  was  not  possible 
otherwise.  They  are 

(1)  A continuum  formulation  with  appropriate  slip  condition  * 
is  physically  plausible  and  can  indeed  be  useful  for  prediction 
purposes  probably  more  reliable  than  experiments. 

(2)  The  surface  conditions  based  on  totally  diffused  reflec- 
tion end  zero  recovery  of  mean  kinetic  energy  is  the  only  correct 
(me  that  can  provide  the  low  surface  pressure  in  the  range  of 
the  available  experimental  data. 

(3)  The  computational  method  with  a reasonable  accuracy  can 
be  useful  in  fluid  dynamics  research  and  can  help  to  resolve  diffi- 
cult fundamental  problems.  It  is  not  siqply  a tool  for  carrying 
out  repetitive  routine  numerical  work  in  engineering  development. 


There  are  important  omissions,  in  the  above  review,  of  many  interes- 
ting and  significant  results  in  the  development  of  co:<putational  methods 
relevant  to  aerospace  applications.  They  are  omitted  here  to  facilitate 
the  presentation  of  the  major  theme  end  hopefully  with  as  little  digres- 
sion as  possible.  While  computational  stability  remains  a problem,  it 
can  generally  be  overcome  with  some  hard  work.  It  should  not  be  permitted 
to  draw  attention  aw  ay  from  the  need  of  reasonably  accurate  computed 
results.  Stable  end  smooth  computational  results  are  encouraging  but  can 
be  very  deceiving.  From  thd  application  point  of  view,  the  question  of 
accuracy  is  crucial.  Accordingly,  the  approach  described  above  to  securo 
"accurate"  formulation  is  of  fundamental  importance,  crude  as  it  is.  How 
such  crude  criteria  may  be  used  and  incorporated  are  desmonstrated  in  this 
chapter.  Much  development  in  this  direction  is  needed.  Some  fundamental 
aspects  should  be  wderstood  and  practical  methods  developed  to  deal  with 
the  various  situations.  Such  problems  will  not  fade  away  because  of  the 
dramatic  advance  in  cocputor  capabilities.  Indeed,  there  are  serious 
problems  thet  will  be  encountered  in  the  efficient  use  of  the  fourth  gen- 
oration  computers  if  any  meaningful  speed  advantages  are  to  be  reaped. 
Therefore,  a few  words  on  the  prospects  of  the  coming  fourth  generation 
confute r will  serve  to  bring  to  conclusion  the  present  review. 
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6.3  Future  Prospects  with  the  Fourth  Generation  Computers 

It  has  been  a constant  allusion  that  faster  and  bigger  computers 
will  provide  the  solution  to  many  difficulties  associated  with  the  numer- 
ical integration  of  partial  differential  equations.  Such  larger  and 
faster  computers  ere  needed  but  they  do  not  provide  the  brutal  force  to 
resolve  all  the  computational  difficulties  without  conscientious  efforts. 

Certain  aspects  of  the  problems  must  be  understood  as  to  their  funda- 
mentals before  being  satisfactorily  dealt  with,  such  as  the  questions  of 
stability  and  accuracy.  Moreover,  the  development  of  the  computer  hard- 
ware has  reached  the  point  through  miniaturisation  the  order  of  magnitude 
improvement  in  the  speed  of  information  processing  cannot  be  expected  as 
was  in  the  past.  The  fourth  generation  computers  promise  to  bring  about 
large  improvement  in  speed  through  "Parallelism"  which  is  very  much  de- 
pendent on  the  sophistication  of  the  software  and  on  the  nature  of  speci- 
fic problems  to  be  solved.  They  bring  complicated  problems  to  the  users 
as  well  as  to  the  manufacturer  of  the  machines. 

"Parallelism"  is  effected  primarily  in  two  different  ways.  Burrough 
Company's  ILLIAC  IV  speeds  up  the  ari thematic  process  by  using  64  ari the- 
matic units,  receiving  the  sane  instruction  from  a conon  command  module 
to  process  simultaneously  64  sets  of  raw  data.  Thus  ari thematic  results 
can  be  "effectively"  obtained  64  times  faster.  This  is  often  referred 
to  as  a "single  Instruction  Multiple  Processor"  machine  (SIMP).  Control 
Data  Corporation's  STAR  (the  STring  ARay  processor)  employs  the  asser&ly 
line  or  "pipe  lino"  technique  in  which  a string  of  data  is  "continuously" 
fed  into  the  "pipe  line"  to  be  processed  by  a standing  instruction.  In 

1 

this  manner,  the  arithematic  unit  does  not  become  idle  when  the  instruc- 
tions are  being  fetched,  decoded,  and  installed  in  place  to  direct  the 
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computation  and  when  the  newly  computed  data  is  being  sent  out  of  the 
arithematic  unit  or  when  the  raw  data  is  being  brought  into  the  arithe- 
aatic  unit.  This  is  often  referred  to  as  the  pipeline  machine,  both 
the  ILLIAC  IV  and  the  STAR  machines  possess  the  virtual  memory  capacity, 
i.e.,  the  machine  will  manage  automatically  the  datm  stored  in  external 
memory  units  for  extending  the  storage  capacity  of  the  machine.  Texas 
Instrument  Corporation's  ASC  machine  (Advanced  Science  Computer)  incor- 
porates both  the  multi -processor  and  the  pipe-line  concept  but  possesses 
no  virtual  memory  capability.  All  these  machines  are  about  to  be  (or 
already)  delivered  by  the  various  manufacturers  and  are  to  become  oper- 
ational shortly. 

ILLIAC  IV  is  most  efficient  when  the  64  arithematic  processors  can 
be  fully  utilized.  Any  vacant  processors  are  simply  idling,  doing  no 
useful  work,  when  the  operation  is  performed  on  less  than  64  sets  of 
data.  Thus  the  demonstration  of  the  speed  of  ILLIAC  IV  vs  currently 
available  computers  is  often  in  terms  of  the  inversion  of  a 64  x 64 
matrix.  STAR  is  most  efficient  when  a large  number  of  raw  data  (the 
long  string  of  data)  is  to  be  processed  through  the  same  operation  so 
that  the  "filling  tine"  of  the  pipe  becooes  negligible  and  the  machine 
will  provide  a 64  time  increase  of  the  effective  speed  since  each  word 
in  the  STAR  contains  64  bits  of  binary  information.  The  ASC  machine  po- 
ssesses intermediate  behavior.  Each  manufacturer  has  developed  powerful 
and  intricate  softwares  to  implement  and  enhance  the  advantages  of  the 
hardware.  But  all  of  them  are  subject  to  the  inherent  limitations  of 
being  a SIMP  or  a pipe- line  machine. 
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For  an/  one  of  these  machines,  a huge  amount  of  data  must  be  stored, 
arranged,  and  retrieved  from  storage  facilities.  They  must  be  done  effic- 
iently, commensurate  with  the  processing  speed  of  the  machine..  Assume 
that  such  can  be  done  for  the  data  in  core  memory,  directly  accessible  to 
the  central  processing  unit  (CPU) , it  cannot  be  done  for  the  data  stored 
in  external  memories.  The  speed  of  a search  operation  and  of  a data  trans- 
mission through  the  interface  to  the  CPU  is  orders  of  magnitude  slower  than 
the  ari thematic  speed  of  the  machine.  If  the  CPU  asks  for  data  in  the  ex- 
ternal storages  too  frequently,  the  CPU  would  be  doing  little  useful  com- 
putation but  transmitting  the  data  in  and  out  of  the  external  memory  units 
under  its  virtual  memory  operation.  If  the  user  should  prefer  to  deprive 
the  machine  of  its  virtual  memory  capability,  then  the  user-programmer 
must  assume  the  responsibility  of  managing  the  data  across  the  interface. 
There  is  an  alternative  solution  of  this  problem  by  expanding  the  core 
memory  of  the  CPU  of  the  computer  to  match  its  processing  speed.  This  is 
infortwately  a very  expensive  proposition.  There  axe  also  other  problems 
of  data  management  in  the  CPU,  probably  not  so  serious  as  the  one  just 
mentioned.  They  are  more  intimately  related  to  the  specific  characteris- 
tics (hardware  and  software)  of  each  computer.  These  ere  the  problems 
which  the  user  cannot  help  very  much  in  its  eventual  solution.  On  the 
other  hand,  these  machines  present  problems  to  the  users,  the  solutions 

a 

of  which  the  manufacturer  of  the  machines  cannot  help. 

Currently  available  computers  are  serial  machines  that  process  and 
advance  the  data  at  one  point  after  another.  Simultaneous  solution  of 
unknowns  at  many  points,  as  is  required  by  implicit  algorithms  is  handled 
through  special  procedures  equivalent  to  matrix  inversion.  If  a program 
designed  for  the  serial  machine  should  be  run  on  the  parallel  computers. 
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no  speed  advantage  will  result.  (Indeed  there  will  be  some  loss.)  The 

64  parallel  processors  of  ILLIAC  IV  will  only  have  one  processor  doing 

, / 

useful  work.  The  STring-ARray  processor  of  STAR  will  operate  in  its 
scalar  node  (versus  the  "vector  node"  for  string  array  processing). 

There  is  not,  and  will  not  be,  a software  that  will  translate  an  existing 
serial  progran  into  a reasonably  efficient  "parallel  program"  for  a 
specific  parallel  machine.  Such  a translation  is  not  a natter  of  trans- 
lating one  language  into  another.  It  is  a natter  of  changing  the  logic 
of  solving  a problem.  It  asks  essentially  for  a new  formulation  for  a 
specific  problem  to  exploit  tho  speed  advantage  offered  by  a specific 
machine.  The  user  is  asked  to  start  anew,  for  each  problen  and  for  a 
specific  ccaputer  and  to  pay  considerable  attention  not  only  to  the  formu- 
lation of  a problem  for  solution  but  also  to  the  storage  of  the  data  in 
the  external  memory  to  match  the  demand  of  the  data  according  to  the  formu- 


lation of  the  problen. 

In  writing  such  a progran  for  use  with  a specific  parallel  conputer, 
it  is  not  a simple  natter  to  take  advantage  of  a successful  serial  progran 
used  with  the  current  serial  machines.  It  may  indeed  be  doubtful,  if 


there  nay  be  any  advantage  under  special  circumstances.  Without  further 
elaboration,  it  nay  be  noted,  even  for  simple  problems,  that: 

1.  An  efficient  serial  algorithm  need  not  lead  to  an  effici- 
ent parallel  algoritha  while  an  inefficient  serial  algorithm 

nay  lead  to  an  efficient  parallel  algoritha. 

2.  A serial  algoritha  that  is  apparently  serial  and  was 
constructed  for  use  with  a serial  conputer  nay  possess  a great 
deal  of  hidden  parallelism  which  nay  be  exploited  to  suit  the 


particular  mode  of  operation  of  a specific  parallel  conputer. 
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3.  Parallel  program  may  behave  quite  differently  in  the 
difference  solution  of  a partial  differential  equation  than 
the  corresponding  serial  program.  Ihe  behavior  refers  to  the 
stability  of  the  computation,  the  rate  of -convergence  to  the 
desired  solution,  and  the  accuracy  of  the  solution. 

Ihe  last  one  is  particularly  important.  It  asks  the  user  to  gain 
as  much  as  possible  the  in derstan cling  of  the  various  fundamental  problems 
of  difference  methods  such  as  stability  and  accuracy.  With  a hotter 
understanding,  it  may  be  hopeful  that  the  years  of  tedious  and  painful 
learning  process  through  trial  end  error  in  developing  the  difference 
techniques  of  the  serial  machines  may  not  be  repeated  or  at  least  may  be 
greatly  reduced. 

For  many  important  practical  problems  the  solution  of  the  Navier- 
Stokes  equations  in  three  spatial  dimensions  will  be  required.  Even  for 
the  steady  state  solution  of  such  problems,  the  confutation  for  a reason- 
ably accurate  solution  will  need  the  speed  and  storage  capacity  promised 
by  these  parallel  computers.  The  complicated  boundary  conditions  do  not 
lend  themselves  to  efficient  parallel  treatments  and  interfere  with  the 
efficient  organization  for  the  parallel  computations  of  the  fluid  flow 
problems.  This  is  in  addition  to  the  fundamental  difficulties  noted 
above.  It  is  much  desired  that  what  has  been  learned  from  the  serial 
machines  may  benefit  the  development  of  computational  programs  that  will 
reap  the  promised  speed  advantage  of  the  parallol  computer.  For  this 
purpose,  it  is  especially  important  to  gain  some  fundamental  understanding 
of  such  complicated  computational  difficulties  special  to  fluid  dynamics. 
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Such  inderstcnding  cannot  be  expected  froa  conputer  scientists  who  have 
their  full  share  of  difficulties  associated  with  the  operation  of  the 
parallel  computers  in  goneral.  Those  wishing  to  solve  the  complicated 
flow  problems  with  the  Navier-Stokes  equations  must  learn  how  to  resolve 


such  difficulties  for  themselves.  Ihe  task  ahead  is  formidable.  The 
potential  reward  is  also  immense. 
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12.  IPONIORIN8  MILITARY  ACTIVITY 


The  V.K.l.  notes  are  being  considered  for  Office  of  Naval  Research 

pifclication  as  a book  in  the  Series  of  Code  438 

Lq^ture  Notes  in  Physics,  Springer  Verlag  Arlington,  Virginia  22217 

.TTU^n^i 

Tne  mathematical  foundation  and  the  various  practical  aspects  of  the  numerical 
solution  of  gas  dynamic  equations  are  critically  reviewed  with  enphasis  on  obtaining 
quantitatively  accurate  solutions  for  application  in  various  engineering  and  sciences 
Computational  stability  rate  of  convergence  and  accuracy  (or  error  estimate)  are  dis- 
cussed. The  promises  and  problems  of  the  4th  generation  computers  are  outlined  with- 
in this  perspective.; 

A‘~ 

- Computational  stability  should  not  be  obtained  at  the  sacrifice  of  the  conver- 
gence rate  to  and  the  accuracy  of  the  final  solution.  With  accuracy  in  mind,  the 
explicit  algorithms  are  likely  preferrable  to  the  iiqplicit  ones.  Strict  conserva- 
tion of  the  difference  formulation  is  recommended  and  exemplified  to  avoid  the  accum- 
ulation of  local  truncation  errors  and  to  facilitate  the  estimate  of  the  errors  in  a 
steady  state  solution.  Illustrative  examples  are  given  including  supersonic  flows 
with  shocks. _ 
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